00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037
00038
00039
00040
00041
00042
00043
00044
00045
00046
00047
00048
00049
00050
00051
00052
00053
00054
00055
00056
00057
00058
00059
00060
00061
00062
00063
00064
00065
00066
00067
00068
00069
00070
00071
00072
00073
00074
00075
00076
00077
00078
00079
00080
00081
00082
00083
00084
00085
00086
00087
00088
00089
00090
00091
00092
00093
00094
00095
00096
00097
00098
00099
00100
00101
00102
00103
00104
00105
00106
00107
00108
00109
00110
00111
00112
00113
00114
00115
00116
00117
00118
00119
00120
00121
00122
00123
00124
00125
00126
00127
00128
00129
00130
00131
00132
00133
00134
00135
00136
00137
00138
00139
00140
00141
00142
00143
00144
00145
00146
00147
00148
00149
00150
00151
00152
00153
00154
00155
00156
00157
00158
00159
00160
00161
00162
00163
00164
00165
00166
00167
00168
00169
00170
00171
00172
00173
00174
00175
00176
00177
00178
00179
00180
00181
00182
00183
00184
00185
00186
00187
00188
00189
00190
00191
00192
00193
00194
00195
00196
00197
00198
00199
00200
00201
00202
00203
00204
00205
00206
00207
00208
00209
00210
00211
00212
00213
00214
00215
00216
00217
00218
00219
00220
00221
00222
00223
00224
00225
00226
00227
00228
00229
00230
00231
00232
00233
00234
00235
00236
00237
00238
00239
00240
00241
00242
00243
00244
00245
00246
00247
00248
00249
00250
00251
00252
00253
00254
00255
00256
00257
00258
00259
00260
00261
00262
00263
00264
00265
00266
00267
00268
00269
00270
00271
00272
00273
00274
00275
00276
00277
00278
00279
00280
00281
00282
00283
00284
00285
00286
00287
00288
00289
00290
00291
00292
00293
00294
00295
00296
00297
00298
00299
00300
00301
00302
00303
00304
00305
00306
00307
00308
00309
00310
00311
00312 #include <stdlib.h>
00313 #include <stdio.h>
00314
00315 #include "TROOT.h"
00316 #include "TMinuit.h"
00317 #include "TMath.h"
00318 #include "TError.h"
00319 #include "TPluginManager.h"
00320 #include "TClass.h"
00321 #include "TInterpreter.h"
00322
00323 TMinuit *gMinuit;
00324
00325 const char charal[29] = " .ABCDEFGHIJKLMNOPQRSTUVWXYZ";
00326
00327 ClassImp(TMinuit)
00328
00329
00330 TMinuit::TMinuit(): TNamed("MINUIT","The Minimization package")
00331 {
00332
00333
00334
00335 if (TMinuit::Class()->IsCallingNew() != TClass::kRealNew) {
00336
00337 fCpnam = 0;
00338 fU = 0;
00339 fAlim = 0;
00340 fBlim = 0;
00341 fPstar = 0;
00342 fGin = 0;
00343 fNvarl = 0;
00344 fNiofex = 0;
00345
00346 fNexofi = 0;
00347 fIpfix = 0;
00348 fErp = 0;
00349 fErn = 0;
00350 fWerr = 0;
00351 fGlobcc = 0;
00352 fX = 0;
00353 fXt = 0;
00354 fDirin = 0;
00355 fXs = 0;
00356 fXts = 0;
00357 fDirins = 0;
00358 fGrd = 0;
00359 fG2 = 0;
00360 fGstep = 0;
00361 fDgrd = 0;
00362 fGrds = 0;
00363 fG2s = 0;
00364 fGsteps = 0;
00365 fPstst = 0;
00366 fPbar = 0;
00367 fPrho = 0;
00368 fWord7 = 0;
00369 fVhmat = 0;
00370 fVthmat = 0;
00371 fP = 0;
00372 fXpt = 0;
00373 fYpt = 0;
00374 fChpt = 0;
00375 fCONTgcc = 0;
00376 fCONTw = 0;
00377 fFIXPyy = 0;
00378 fGRADgf = 0;
00379 fHESSyy = 0;
00380 fIMPRdsav = 0;
00381 fIMPRy = 0;
00382 fMATUvline = 0;
00383 fMIGRflnu = 0;
00384 fMIGRstep = 0;
00385 fMIGRgs = 0;
00386 fMIGRvg = 0;
00387 fMIGRxxs = 0;
00388 fMNOTxdev = 0;
00389 fMNOTw = 0;
00390 fMNOTgcc = 0;
00391 fPSDFs = 0;
00392 fSEEKxmid = 0;
00393 fSEEKxbest = 0;
00394 fSIMPy = 0;
00395 fVERTq = 0;
00396 fVERTs = 0;
00397 fVERTpp = 0;
00398 fCOMDplist = 0;
00399 fPARSplist = 0;
00400
00401 fUp = 0;
00402 fEpsi = 0;
00403 fApsi = 0;
00404 fXmidcr = 0;
00405 fYmidcr = 0;
00406 fXdircr = 0;
00407 fYdircr = 0;
00408
00409 fStatus = 0;
00410 fEmpty = 0;
00411 fObjectFit = 0;
00412 fMethodCall = 0;
00413 fPlot = 0;
00414 fGraphicsMode = kTRUE;
00415
00416 } else {
00417 BuildArrays(25);
00418
00419 fUp = 0;
00420 fEpsi = 0;
00421 fApsi = 0;
00422 fXmidcr = 0;
00423 fYmidcr = 0;
00424 fXdircr = 0;
00425 fYdircr = 0;
00426
00427 fStatus = 0;
00428 fEmpty = 0;
00429 fObjectFit = 0;
00430 fMethodCall = 0;
00431 fPlot = 0;
00432 fGraphicsMode = kTRUE;
00433 SetMaxIterations();
00434 mninit(5,6,7);
00435 }
00436
00437 fFCN = 0;
00438 gMinuit = this;
00439 gROOT->GetListOfSpecials()->Add(gMinuit);
00440
00441 }
00442
00443
00444 TMinuit::TMinuit(Int_t maxpar): TNamed("MINUIT","The Minimization package")
00445 {
00446
00447
00448
00449
00450
00451 fFCN = 0;
00452
00453 BuildArrays(maxpar);
00454
00455 fStatus = 0;
00456 fEmpty = 0;
00457 fObjectFit = 0;
00458 fMethodCall = 0;
00459 fPlot = 0;
00460 fGraphicsMode = kTRUE;
00461 SetMaxIterations();
00462
00463 mninit(5,6,7);
00464 gMinuit = this;
00465 gROOT->GetListOfSpecials()->Add(gMinuit);
00466 }
00467
00468
00469 TMinuit::TMinuit(const TMinuit &minuit) : TNamed(minuit)
00470 {
00471
00472
00473 Error("TMinuit", "can not copy construct TMinuit");
00474 }
00475
00476
00477 TMinuit::~TMinuit()
00478 {
00479
00480
00481
00482 DeleteArrays();
00483 delete fPlot;
00484 delete fMethodCall;
00485 gROOT->GetListOfSpecials()->Remove(this);
00486 if (gMinuit == this) gMinuit = 0;
00487 }
00488
00489
00490 void TMinuit::BuildArrays(Int_t maxpar)
00491 {
00492
00493
00494
00495 fMaxpar = 25;
00496 if (maxpar >= fMaxpar) fMaxpar = maxpar+1;
00497 fMaxpar1= fMaxpar*(fMaxpar+1);
00498 fMaxpar2= 2*fMaxpar;
00499 fMaxpar5= fMaxpar1/2;
00500 fMaxcpt = 101;
00501 fCpnam = new TString[fMaxpar2];
00502 fU = new Double_t[fMaxpar2];
00503 fAlim = new Double_t[fMaxpar2];
00504 fBlim = new Double_t[fMaxpar2];
00505 fPstar = new Double_t[fMaxpar2];
00506 fGin = new Double_t[fMaxpar2];
00507 fNvarl = new Int_t[fMaxpar2];
00508 fNiofex = new Int_t[fMaxpar2];
00509
00510 fNexofi = new Int_t[fMaxpar];
00511 fIpfix = new Int_t[fMaxpar];
00512 fErp = new Double_t[fMaxpar];
00513 fErn = new Double_t[fMaxpar];
00514 fWerr = new Double_t[fMaxpar];
00515 fGlobcc = new Double_t[fMaxpar];
00516 fX = new Double_t[fMaxpar];
00517 fXt = new Double_t[fMaxpar];
00518 fDirin = new Double_t[fMaxpar];
00519 fXs = new Double_t[fMaxpar];
00520 fXts = new Double_t[fMaxpar];
00521 fDirins = new Double_t[fMaxpar];
00522 fGrd = new Double_t[fMaxpar];
00523 fG2 = new Double_t[fMaxpar];
00524 fGstep = new Double_t[fMaxpar];
00525 fDgrd = new Double_t[fMaxpar];
00526 fGrds = new Double_t[fMaxpar];
00527 fG2s = new Double_t[fMaxpar];
00528 fGsteps = new Double_t[fMaxpar];
00529 fPstst = new Double_t[fMaxpar];
00530 fPbar = new Double_t[fMaxpar];
00531 fPrho = new Double_t[fMaxpar];
00532 fWord7 = new Double_t[fMaxpar];
00533 fVhmat = new Double_t[fMaxpar5];
00534 fVthmat = new Double_t[fMaxpar5];
00535 fP = new Double_t[fMaxpar1];
00536 fXpt = new Double_t[fMaxcpt];
00537 fYpt = new Double_t[fMaxcpt];
00538 fChpt = new char[fMaxcpt+1];
00539
00540
00541 fCONTgcc = new Double_t[fMaxpar];
00542 fCONTw = new Double_t[fMaxpar];
00543 fFIXPyy = new Double_t[fMaxpar];
00544 fGRADgf = new Double_t[fMaxpar];
00545 fHESSyy = new Double_t[fMaxpar];
00546 fIMPRdsav = new Double_t[fMaxpar];
00547 fIMPRy = new Double_t[fMaxpar];
00548 fMATUvline = new Double_t[fMaxpar];
00549 fMIGRflnu = new Double_t[fMaxpar];
00550 fMIGRstep = new Double_t[fMaxpar];
00551 fMIGRgs = new Double_t[fMaxpar];
00552 fMIGRvg = new Double_t[fMaxpar];
00553 fMIGRxxs = new Double_t[fMaxpar];
00554 fMNOTxdev = new Double_t[fMaxpar];
00555 fMNOTw = new Double_t[fMaxpar];
00556 fMNOTgcc = new Double_t[fMaxpar];
00557 fPSDFs = new Double_t[fMaxpar];
00558 fSEEKxmid = new Double_t[fMaxpar];
00559 fSEEKxbest = new Double_t[fMaxpar];
00560 fSIMPy = new Double_t[fMaxpar];
00561 fVERTq = new Double_t[fMaxpar];
00562 fVERTs = new Double_t[fMaxpar];
00563 fVERTpp = new Double_t[fMaxpar];
00564 fCOMDplist = new Double_t[fMaxpar];
00565 fPARSplist = new Double_t[fMaxpar];
00566
00567 for (int i = 0; i < fMaxpar; i++) {
00568 fErp[i] = 0;
00569 fErn[i] = 0;
00570 }
00571 }
00572
00573
00574
00575 TObject *TMinuit::Clone(const char *newname) const
00576 {
00577
00578
00579 TMinuit *named = (TMinuit*)TNamed::Clone(newname);
00580 named->fFCN=fFCN;
00581 return named;
00582 }
00583
00584
00585
00586 Int_t TMinuit::Command(const char *command)
00587 {
00588
00589
00590
00591
00592
00593
00594
00595
00596
00597
00598
00599
00600
00601
00602
00603
00604
00605
00606
00607
00608
00609
00610
00611
00612 Int_t status = 0;
00613 mncomd(command,status);
00614 return status;
00615 }
00616
00617
00618 TObject *TMinuit::Contour(Int_t npoints, Int_t pa1, Int_t pa2)
00619 {
00620
00621
00622
00623
00624
00625
00626
00627
00628
00629
00630
00631
00632
00633
00634
00635
00636
00637
00638
00639 if (npoints<4) {
00640
00641 fStatus= 2;
00642 return (TObject *)0;
00643 }
00644 Int_t npfound;
00645 Double_t *xcoor = new Double_t[npoints+1];
00646 Double_t *ycoor = new Double_t[npoints+1];
00647 mncont(pa1,pa2,npoints,xcoor,ycoor,npfound);
00648 if (npfound<4) {
00649
00650 Warning("Contour","Cannot find more than 4 points, no TGraph returned");
00651 fStatus= (npfound==0 ? 1 : npfound);
00652 delete [] xcoor;
00653 delete [] ycoor;
00654 return (TObject *)0;
00655 }
00656 if (npfound!=npoints) {
00657
00658 Warning("Contour","Returning a TGraph with %d points only",npfound);
00659 npoints = npfound;
00660 }
00661 fStatus=0;
00662
00663 xcoor[npoints] = xcoor[0];
00664 ycoor[npoints] = ycoor[0];
00665 TObject *gr = 0;
00666 TPluginHandler *h;
00667 if ((h = gROOT->GetPluginManager()->FindHandler("TMinuitGraph"))) {
00668 if (h->LoadPlugin() != -1)
00669 gr = (TObject*)h->ExecPlugin(3,npoints+1,xcoor,ycoor);
00670 }
00671 delete [] xcoor;
00672 delete [] ycoor;
00673 return gr;
00674 }
00675
00676
00677 Int_t TMinuit::DefineParameter( Int_t parNo, const char *name, Double_t initVal, Double_t initErr, Double_t lowerLimit, Double_t upperLimit )
00678 {
00679
00680
00681 Int_t err;
00682
00683 TString sname = name;
00684 mnparm( parNo, sname, initVal, initErr, lowerLimit, upperLimit, err);
00685
00686 return err;
00687 }
00688
00689
00690 void TMinuit::DeleteArrays()
00691 {
00692
00693
00694 if (fEmpty) return;
00695 delete [] fCpnam;
00696 delete [] fU;
00697 delete [] fAlim;
00698 delete [] fBlim;
00699 delete [] fErp;
00700 delete [] fErn;
00701 delete [] fWerr;
00702 delete [] fGlobcc;
00703 delete [] fNvarl;
00704 delete [] fNiofex;
00705 delete [] fNexofi;
00706 delete [] fX;
00707 delete [] fXt;
00708 delete [] fDirin;
00709 delete [] fXs;
00710 delete [] fXts;
00711 delete [] fDirins;
00712 delete [] fGrd;
00713 delete [] fG2;
00714 delete [] fGstep;
00715 delete [] fGin;
00716 delete [] fDgrd;
00717 delete [] fGrds;
00718 delete [] fG2s;
00719 delete [] fGsteps;
00720 delete [] fIpfix;
00721 delete [] fVhmat;
00722 delete [] fVthmat;
00723 delete [] fP;
00724 delete [] fPstar;
00725 delete [] fPstst;
00726 delete [] fPbar;
00727 delete [] fPrho;
00728 delete [] fWord7;
00729 delete [] fXpt;
00730 delete [] fYpt;
00731 delete [] fChpt;
00732
00733 delete [] fCONTgcc;
00734 delete [] fCONTw;
00735 delete [] fFIXPyy;
00736 delete [] fGRADgf;
00737 delete [] fHESSyy;
00738 delete [] fIMPRdsav;
00739 delete [] fIMPRy;
00740 delete [] fMATUvline;
00741 delete [] fMIGRflnu;
00742 delete [] fMIGRstep;
00743 delete [] fMIGRgs;
00744 delete [] fMIGRvg;
00745 delete [] fMIGRxxs;
00746 delete [] fMNOTxdev;
00747 delete [] fMNOTw;
00748 delete [] fMNOTgcc;
00749 delete [] fPSDFs;
00750 delete [] fSEEKxmid;
00751 delete [] fSEEKxbest;
00752 delete [] fSIMPy;
00753 delete [] fVERTq;
00754 delete [] fVERTs;
00755 delete [] fVERTpp;
00756 delete [] fCOMDplist;
00757 delete [] fPARSplist;
00758
00759 fEmpty = 1;
00760 }
00761
00762
00763 Int_t TMinuit::Eval(Int_t npar, Double_t *grad, Double_t &fval, Double_t *par, Int_t flag)
00764 {
00765
00766
00767
00768
00769
00770
00771
00772
00773
00774
00775
00776
00777
00778
00779
00780
00781
00782
00783
00784
00785
00786
00787
00788
00789
00790
00791
00792
00793
00794
00795
00796
00797
00798
00799
00800
00801
00802
00803 if (fFCN) (*fFCN)(npar,grad,fval,par,flag);
00804 return 0;
00805 }
00806
00807
00808 Int_t TMinuit::FixParameter( Int_t parNo)
00809 {
00810
00811
00812 Int_t err;
00813 Double_t tmp[1];
00814 tmp[0] = parNo+1;
00815
00816 mnexcm( "FIX", tmp, 1, err );
00817
00818 return err;
00819 }
00820
00821
00822 Int_t TMinuit::GetParameter( Int_t parNo, Double_t ¤tValue, Double_t ¤tError ) const
00823 {
00824
00825 Int_t err;
00826 TString name;
00827 Double_t bnd1, bnd2;
00828
00829 mnpout( parNo, name, currentValue, currentError, bnd1, bnd2, err );
00830
00831 return err;
00832 }
00833
00834
00835 Int_t TMinuit::GetNumFixedPars() const
00836 {
00837
00838
00839 return fNpfix;
00840 }
00841
00842
00843 Int_t TMinuit::GetNumFreePars() const
00844 {
00845
00846
00847 return fNpar;
00848 }
00849
00850
00851 Int_t TMinuit::GetNumPars() const
00852 {
00853
00854
00855
00856 return fNpar + fNpfix;
00857 }
00858
00859
00860 Int_t TMinuit::Migrad()
00861 {
00862
00863 Int_t err;
00864 Double_t tmp[1];
00865 tmp[0] = 0;
00866
00867 mnexcm( "MIGRAD", tmp, 0, err );
00868
00869 return err;
00870 }
00871
00872
00873 Int_t TMinuit::Release( Int_t parNo)
00874 {
00875
00876
00877 Int_t err;
00878 Double_t tmp[1];
00879 tmp[0] = parNo+1;
00880
00881 mnexcm( "RELEASE", tmp, 1, err );
00882
00883 return err;
00884 }
00885
00886
00887 Int_t TMinuit::SetErrorDef( Double_t up )
00888 {
00889
00890
00891 Int_t err;
00892
00893 mnexcm( "SET ERRDEF", &up, 1, err );
00894
00895 return err;
00896 }
00897
00898
00899 void TMinuit::SetFCN(void (*fcn)(Int_t &, Double_t *, Double_t &f, Double_t *, Int_t))
00900 {
00901
00902
00903
00904 fFCN = fcn;
00905 }
00906
00907
00908 void InteractiveFCNm(Int_t &npar, Double_t *gin, Double_t &f, Double_t *u, Int_t flag)
00909 {
00910
00911
00912
00913 TMethodCall *m = gMinuit->GetMethodCall();
00914 if (!m) return;
00915
00916 Long_t args[5];
00917 args[0] = (Long_t)∦
00918 args[1] = (Long_t)gin;
00919 args[2] = (Long_t)&f;
00920 args[3] = (Long_t)u;
00921 args[4] = (Long_t)flag;
00922 m->SetParamPtrs(args);
00923 Double_t result;
00924 m->Execute(result);
00925 }
00926
00927
00928 void TMinuit::SetFCN(void *fcn)
00929 {
00930
00931
00932
00933
00934
00935 if (!fcn) return;
00936
00937 const char *funcname = gCint->Getp2f2funcname(fcn);
00938 if (funcname) {
00939 fMethodCall = new TMethodCall();
00940 fMethodCall->InitWithPrototype(funcname,"Int_t&,Double_t*,Double_t&,Double_t*,Int_t");
00941 }
00942 fFCN = InteractiveFCNm;
00943 gMinuit = this;
00944 }
00945
00946
00947 Int_t TMinuit::SetPrintLevel( Int_t printLevel )
00948 {
00949
00950
00951
00952
00953 Int_t err;
00954 Double_t tmp[1];
00955 tmp[0] = printLevel;
00956
00957 mnexcm( "SET PRINT", tmp, 1, err );
00958
00959 if (printLevel <=-1) mnexcm("SET NOWarnings",tmp,0,err);
00960
00961 return err;
00962 }
00963
00964
00965 void TMinuit::mnamin()
00966 {
00967
00968
00969
00970
00971
00972
00973
00974
00975 Double_t fnew;
00976 Int_t nparx;
00977
00978 nparx = fNpar;
00979 if (fISW[4] >= 1) {
00980 Printf(" FIRST CALL TO USER FUNCTION AT NEW START POINT, WITH IFLAG=4.");
00981 }
00982 mnexin(fX);
00983 Eval(nparx, fGin, fnew, fU, 4); ++fNfcn;
00984 fAmin = fnew;
00985 fEDM = fBigedm;
00986 }
00987
00988
00989 void TMinuit::mnbins(Double_t a1, Double_t a2, Int_t naa, Double_t &bl, Double_t &bh, Int_t &nb, Double_t &bwid)
00990 {
00991
00992
00993
00994
00995
00996
00997
00998
00999
01000
01001 Double_t awid,ah, al, sigfig, sigrnd, alb;
01002 Int_t kwid, lwid, na=0, log_;
01003
01004 al = TMath::Min(a1,a2);
01005 ah = TMath::Max(a1,a2);
01006 if (al == ah) ah = al + 1;
01007
01008
01009 if (naa == -1) goto L150;
01010 L10:
01011 na = naa - 1;
01012 if (na < 1) na = 1;
01013
01014
01015 L20:
01016 awid = (ah-al) / Double_t(na);
01017 log_ = Int_t(TMath::Log10(awid));
01018 if (awid <= 1) --log_;
01019 sigfig = awid*TMath::Power(10, -log_);
01020
01021 if (sigfig > 2) goto L40;
01022 sigrnd = 2;
01023 goto L100;
01024 L40:
01025 if (sigfig > 2.5) goto L50;
01026 sigrnd = 2.5;
01027 goto L100;
01028 L50:
01029 if (sigfig > 5) goto L60;
01030 sigrnd = 5;
01031 goto L100;
01032 L60:
01033 sigrnd = 1;
01034 ++log_;
01035 L100:
01036 bwid = sigrnd*TMath::Power(10, log_);
01037 goto L200;
01038
01039 L150:
01040 if (bwid <= 0) goto L10;
01041 L200:
01042 alb = al / bwid;
01043 lwid = Int_t(alb);
01044 if (alb < 0) --lwid;
01045 bl = bwid*Double_t(lwid);
01046 alb = ah / bwid + 1;
01047 kwid = Int_t(alb);
01048 if (alb < 0) --kwid;
01049 bh = bwid*Double_t(kwid);
01050 nb = kwid - lwid;
01051 if (naa > 5) goto L240;
01052 if (naa == -1) return;
01053
01054 if (naa > 1 || nb == 1) return;
01055 bwid *= 2;
01056 nb = 1;
01057 return;
01058 L240:
01059 if (nb << 1 != naa) return;
01060 ++na;
01061 goto L20;
01062 }
01063
01064
01065 void TMinuit::mncalf(Double_t *pvec, Double_t &ycalf)
01066 {
01067
01068
01069
01070
01071
01072
01073
01074
01075 Int_t ndex, i, j, m, n, nparx;
01076 Double_t denom, f;
01077
01078 nparx = fNpar;
01079 mninex(&pvec[0]);
01080 Eval(nparx, fGin, f, fU, 4); ++fNfcn;
01081 for (i = 1; i <= fNpar; ++i) {
01082 fGrd[i-1] = 0;
01083 for (j = 1; j <= fNpar; ++j) {
01084 m = TMath::Max(i,j);
01085 n = TMath::Min(i,j);
01086 ndex = m*(m-1) / 2 + n;
01087 fGrd[i-1] += fVthmat[ndex-1]*(fXt[j-1] - pvec[j-1]);
01088 }
01089 }
01090 denom = 0;
01091 for (i = 1; i <= fNpar; ++i) {denom += fGrd[i-1]*(fXt[i-1] - pvec[i-1]); }
01092 if (denom <= 0) {
01093 fDcovar = 1;
01094 fISW[1] = 0;
01095 denom = 1;
01096 }
01097 ycalf = (f - fApsi) / denom;
01098 }
01099
01100
01101 void TMinuit::mncler()
01102 {
01103
01104
01105
01106
01107
01108 Int_t i;
01109
01110 fNpfix = 0;
01111 fNu = 0;
01112 fNpar = 0;
01113 fNfcn = 0;
01114 fNwrmes[0] = 0;
01115 fNwrmes[1] = 0;
01116 for (i = 1; i <= fMaxext; ++i) {
01117 fU[i-1] = 0;
01118 fCpnam[i-1] = fCundef;
01119 fNvarl[i-1] = -1;
01120 fNiofex[i-1] = 0;
01121 }
01122 mnrset(1);
01123 fCfrom = "CLEAR ";
01124 fNfcnfr = fNfcn;
01125 fCstatu = "UNDEFINED ";
01126 fLnolim = kTRUE;
01127 fLphead = kTRUE;
01128 }
01129
01130
01131 void TMinuit::mncntr(Int_t ike1, Int_t ike2, Int_t &ierrf)
01132 {
01133
01134
01135
01136
01137
01138
01139 static TString clabel = "0123456789ABCDEFGHIJ";
01140
01141
01142 Double_t d__1, d__2;
01143 Double_t fcna[115], fcnb[115], contur[20];
01144 Double_t ylabel, fmn, fmx, xlo, ylo, xup, yup;
01145 Double_t devs, xsav, ysav, bwidx, bwidy, unext, ff, xb4;
01146 Int_t i, ngrid, ixmid, nparx, ix, nx, ny, ki1, ki2, ixzero, iy, ics;
01147 TString chmid, chln, chzero;
01148
01149 Int_t ke1 = ike1+1;
01150 Int_t ke2 = ike2+1;
01151 if (ke1 <= 0 || ke2 <= 0) goto L1350;
01152 if (ke1 > fNu || ke2 > fNu) goto L1350;
01153 ki1 = fNiofex[ke1-1];
01154 ki2 = fNiofex[ke2-1];
01155 if (ki1 <= 0 || ki2 <= 0) goto L1350;
01156 if (ki1 == ki2) goto L1350;
01157
01158 if (fISW[1] < 1) {
01159 mnhess();
01160 mnwerr();
01161 }
01162 nparx = fNpar;
01163 xsav = fU[ke1-1];
01164 ysav = fU[ke2-1];
01165 devs = fWord7[2];
01166 if (devs <= 0) devs = 2;
01167 xlo = fU[ke1-1] - devs*fWerr[ki1-1];
01168 xup = fU[ke1-1] + devs*fWerr[ki1-1];
01169 ylo = fU[ke2-1] - devs*fWerr[ki2-1];
01170 yup = fU[ke2-1] + devs*fWerr[ki2-1];
01171 ngrid = Int_t(fWord7[3]);
01172 if (ngrid <= 0) {
01173 ngrid = 25;
01174
01175 nx = TMath::Min(fNpagwd - 15,ngrid);
01176
01177 ny = TMath::Min(fNpagln - 7,ngrid);
01178 } else {
01179 nx = ngrid;
01180 ny = ngrid;
01181 }
01182 if (nx < 11) nx = 11;
01183 if (ny < 11) ny = 11;
01184 if (nx >= 115) nx = 114;
01185
01186
01187 if (fNvarl[ke1-1] > 1) {
01188 if (xlo < fAlim[ke1-1]) xlo = fAlim[ke1-1];
01189 if (xup > fBlim[ke1-1]) xup = fBlim[ke1-1];
01190 }
01191 if (fNvarl[ke2-1] > 1) {
01192 if (ylo < fAlim[ke2-1]) ylo = fAlim[ke2-1];
01193 if (yup > fBlim[ke2-1]) yup = fBlim[ke2-1];
01194 }
01195 bwidx = (xup - xlo) / Double_t(nx);
01196 bwidy = (yup - ylo) / Double_t(ny);
01197 ixmid = Int_t(((xsav - xlo)*Double_t(nx) / (xup - xlo)) + 1);
01198 if (fAmin == fUndefi) mnamin();
01199
01200 for (i = 1; i <= 20; ++i) { contur[i-1] = fAmin + fUp*(i-1)*(i-1); }
01201 contur[0] += fUp*.01;
01202
01203 fU[ke2-1] = yup;
01204 ixzero = 0;
01205 xb4 = 1;
01206
01207 chmid.Resize(nx+1);
01208 chzero.Resize(nx+1);
01209 chln.Resize(nx+1);
01210 for (ix = 1; ix <= nx + 1; ++ix) {
01211 fU[ke1-1] = xlo + Double_t(ix-1)*bwidx;
01212 Eval(nparx, fGin, ff, fU, 4);
01213 fcnb[ix-1] = ff;
01214 if (xb4 < 0 && fU[ke1-1] > 0) ixzero = ix - 1;
01215 xb4 = fU[ke1-1];
01216 chmid[ix-1] = '*';
01217 chzero[ix-1] = '-';
01218 }
01219 Printf(" Y-AXIS: PARAMETER %3d: %s",ke2,(const char*)fCpnam[ke2-1]);
01220 if (ixzero > 0) {
01221 chzero[ixzero-1] = '+';
01222 chln = " ";
01223 Printf(" X=0");
01224 }
01225
01226 for (iy = 1; iy <= ny; ++iy) {
01227 unext = fU[ke2-1] - bwidy;
01228
01229 chln = " ";
01230
01231 chln.Resize(nx+1);
01232 chln[ixmid-1] = '*';
01233 if (ixzero != 0) chln[ixzero-1] = ':';
01234 if (fU[ke2-1] > ysav && unext < ysav) chln = chmid;
01235 if (fU[ke2-1] > 0 && unext < 0) chln = chzero;
01236 fU[ke2-1] = unext;
01237 ylabel = fU[ke2-1] + bwidy*.5;
01238
01239 for (ix = 1; ix <= nx + 1; ++ix) {
01240 fcna[ix-1] = fcnb[ix-1];
01241 fU[ke1-1] = xlo + Double_t(ix-1)*bwidx;
01242 Eval(nparx, fGin, ff, fU, 4);
01243 fcnb[ix-1] = ff;
01244 }
01245
01246 for (ix = 1; ix <= nx; ++ix) {
01247 d__1 = TMath::Max(fcna[ix-1],fcnb[ix-1]),
01248 d__2 = TMath::Max(fcna[ix],fcnb[ix]);
01249 fmx = TMath::Max(d__1,d__2);
01250 d__1 = TMath::Min(fcna[ix-1],fcnb[ix-1]),
01251 d__2 = TMath::Min(fcna[ix],fcnb[ix]);
01252 fmn = TMath::Min(d__1,d__2);
01253 for (ics = 1; ics <= 20; ++ics) {
01254 if (contur[ics-1] > fmn) goto L240;
01255 }
01256 continue;
01257 L240:
01258 if (contur[ics-1] < fmx) chln[ix-1] = clabel[ics-1];
01259 }
01260
01261 Printf(" %12.4g %s",ylabel,(const char*)chln);
01262 }
01263
01264 chln = " ";
01265 chln(0,1) = 'I';
01266 chln(ixmid-1,1) = 'I';
01267 chln(nx-1,1) = 'I';
01268 Printf(" %s",(const char*)chln);
01269
01270
01271 chln = " ";
01272 if (nx <= 26) {
01273 Printf(" %12.4g%s%12.4g",xlo,(const char*)chln,xup);
01274 Printf(" %s%12.4g",(const char*)chln,xsav);
01275 } else {
01276 Printf(" %12.4g%s%12.4g%s%12.4g",xlo,(const char*)chln,xsav,(const char*)chln,xup);
01277 }
01278 Printf(" X-AXIS: PARAMETER %3d %s ONE COLUMN=%12.4g"
01279 ,ke1,(const char*)fCpnam[ke1-1],bwidx);
01280 Printf(" FUNCTION VALUES: F(I)=%12.4g +%12.4g *I**2",fAmin,fUp);
01281
01282 fU[ke1-1] = xsav;
01283 fU[ke2-1] = ysav;
01284 ierrf = 0;
01285 return;
01286 L1350:
01287 Printf(" INVALID PARAMETER NUMBER(S) REQUESTED. IGNORED.");
01288 ierrf = 1;
01289 }
01290
01291
01292 void TMinuit::mncomd(const char *crdbin, Int_t &icondn)
01293 {
01294
01295
01296
01297
01298
01299
01300
01301
01302
01303
01304
01305
01306
01307
01308
01309
01310
01311
01312
01313
01314
01315
01316
01317 Int_t ierr, ipos, i, llist, lenbuf, lnc;
01318 Bool_t leader;
01319 TString comand, crdbuf, ctemp;
01320
01321 crdbuf = crdbin;
01322 crdbuf.ToUpper();
01323 lenbuf = crdbuf.Length();
01324 icondn = 0;
01325
01326 leader = kTRUE;
01327 ipos = 1;
01328 for (i = 1; i <= TMath::Min(20,lenbuf); ++i) {
01329 if (crdbuf[i-1] == '\'') break;
01330 if (crdbuf[i-1] == ' ') {
01331 if (leader) ++ipos;
01332 continue;
01333 }
01334 leader = kFALSE;
01335 }
01336
01337
01338 if (ipos > lenbuf) {
01339 Printf(" BLANK COMMAND IGNORED.");
01340 icondn = 1;
01341 return;
01342 }
01343
01344
01345 if (crdbuf(ipos-1,3) == "PAR") {
01346 icondn = 5;
01347 fLphead = kTRUE;
01348 return;
01349 }
01350
01351 if (crdbuf(ipos-1,3) == "SET INP") {
01352 icondn = 6;
01353 fLphead = kTRUE;
01354 return;
01355 }
01356
01357 if (crdbuf(ipos-1,7) == "SET TIT") {
01358 icondn = 7;
01359 fLphead = kTRUE;
01360 return;
01361 }
01362
01363 if (crdbuf(ipos-1,7) == "SET COV") {
01364 icondn = 8;
01365 fLphead = kTRUE;
01366 return;
01367 }
01368
01369 ctemp = crdbuf(ipos-1,lenbuf-ipos+1);
01370 mncrck(ctemp, 20, comand, lnc, fMaxpar, fCOMDplist, llist, ierr, fIsyswr);
01371 if (ierr > 0) {
01372 Printf(" COMMAND CANNOT BE INTERPRETED");
01373 icondn = 2;
01374 return;
01375 }
01376
01377 mnexcm(comand.Data(), fCOMDplist, llist, ierr);
01378 icondn = ierr;
01379 }
01380
01381
01382 void TMinuit::mncont(Int_t ike1, Int_t ike2, Int_t nptu, Double_t *xptu, Double_t *yptu, Int_t &ierrf)
01383 {
01384
01385
01386
01387
01388
01389
01390
01391
01392
01393
01394
01395
01396
01397
01398
01399 Int_t i__1;
01400
01401
01402 Double_t d__1, d__2;
01403 Double_t dist, xdir, ydir, aopt, u1min, u2min;
01404 Double_t abest, scalx, scaly;
01405 Double_t a1, a2, val2mi, val2pl, dc, sclfac, bigdis, sigsav;
01406 Int_t nall, iold, line, mpar, ierr, inew, move, next, i, j, nfcol, iercr;
01407 Int_t idist=0, npcol, kints, i2, i1, lr, nfcnco=0, ki1, ki2, ki3, ke3;
01408 Int_t nowpts, istrav, nfmxin, isw2, isw4;
01409 Bool_t ldebug;
01410
01411
01412 Int_t ke1 = ike1+1;
01413 Int_t ke2 = ike2+1;
01414 ldebug = fIdbg[6] >= 1;
01415 if (ke1 <= 0 || ke2 <= 0) goto L1350;
01416 if (ke1 > fNu || ke2 > fNu) goto L1350;
01417 ki1 = fNiofex[ke1-1];
01418 ki2 = fNiofex[ke2-1];
01419 if (ki1 <= 0 || ki2 <= 0) goto L1350;
01420 if (ki1 == ki2) goto L1350;
01421 if (nptu < 4) goto L1400;
01422
01423 nfcnco = fNfcn;
01424 fNfcnmx = (nptu + 5)*100*(fNpar + 1);
01425
01426 mncuve();
01427 u1min = fU[ke1-1];
01428 u2min = fU[ke2-1];
01429 ierrf = 0;
01430 fCfrom = "MNContour ";
01431 fNfcnfr = nfcnco;
01432 if (fISW[4] >= 0) {
01433 Printf(" START MNCONTOUR CALCULATION OF %4d POINTS ON CONTOUR.",nptu);
01434 if (fNpar > 2) {
01435 if (fNpar == 3) {
01436 ki3 = 6 - ki1 - ki2;
01437 ke3 = fNexofi[ki3-1];
01438 Printf(" EACH POINT IS A MINIMUM WITH RESPECT TO PARAMETER %3d %s",ke3,(const char*)fCpnam[ke3-1]);
01439 } else {
01440 Printf(" EACH POINT IS A MINIMUM WITH RESPECT TO THE OTHER %3d VARIABLE PARAMETERS.",fNpar - 2);
01441 }
01442 }
01443 }
01444
01445
01446
01447 mnmnot(ke1, ke2, val2pl, val2mi);
01448 if (fErn[ki1-1] == fUndefi) {
01449 xptu[0] = fAlim[ke1-1];
01450 mnwarn("W", "MNContour ", "Contour squeezed by parameter limits.");
01451 } else {
01452 if (fErn[ki1-1] >= 0) goto L1500;
01453 xptu[0] = u1min + fErn[ki1-1];
01454 }
01455 yptu[0] = val2mi;
01456
01457 if (fErp[ki1-1] == fUndefi) {
01458 xptu[2] = fBlim[ke1-1];
01459 mnwarn("W", "MNContour ", "Contour squeezed by parameter limits.");
01460 } else {
01461 if (fErp[ki1-1] <= 0) goto L1500;
01462 xptu[2] = u1min + fErp[ki1-1];
01463 }
01464 yptu[2] = val2pl;
01465 scalx = 1 / (xptu[2] - xptu[0]);
01466
01467 mnmnot(ke2, ke1, val2pl, val2mi);
01468 if (fErn[ki2-1] == fUndefi) {
01469 yptu[1] = fAlim[ke2-1];
01470 mnwarn("W", "MNContour ", "Contour squeezed by parameter limits.");
01471 } else {
01472 if (fErn[ki2-1] >= 0) goto L1500;
01473 yptu[1] = u2min + fErn[ki2-1];
01474 }
01475 xptu[1] = val2mi;
01476 if (fErp[ki2-1] == fUndefi) {
01477 yptu[3] = fBlim[ke2-1];
01478 mnwarn("W", "MNContour ", "Contour squeezed by parameter limits.");
01479 } else {
01480 if (fErp[ki2-1] <= 0) goto L1500;
01481 yptu[3] = u2min + fErp[ki2-1];
01482 }
01483 xptu[3] = val2pl;
01484 scaly = 1 / (yptu[3] - yptu[1]);
01485 nowpts = 4;
01486 next = 5;
01487 if (ldebug) {
01488 Printf(" Plot of four points found by MINOS");
01489 fXpt[0] = u1min;
01490 fYpt[0] = u2min;
01491 fChpt[0] = ' ';
01492
01493 nall = TMath::Min(nowpts + 1,101);
01494 for (i = 2; i <= nall; ++i) {
01495 fXpt[i-1] = xptu[i-2];
01496 fYpt[i-1] = yptu[i-2];
01497 }
01498 sprintf(fChpt,"%s"," ABCD");
01499 mnplot(fXpt, fYpt, fChpt, nall, fNpagwd, fNpagln);
01500 }
01501
01502
01503 isw2 = fISW[1];
01504 isw4 = fISW[3];
01505 sigsav = fEDM;
01506 istrav = fIstrat;
01507 dc = fDcovar;
01508 fApsi = fEpsi*.5;
01509 abest = fAmin;
01510 mpar = fNpar;
01511 nfmxin = fNfcnmx;
01512 for (i = 1; i <= mpar; ++i) { fXt[i-1] = fX[i-1]; }
01513 i__1 = mpar*(mpar + 1) / 2;
01514 for (j = 1; j <= i__1; ++j) { fVthmat[j-1] = fVhmat[j-1]; }
01515 for (i = 1; i <= mpar; ++i) {
01516 fCONTgcc[i-1] = fGlobcc[i-1];
01517 fCONTw[i-1] = fWerr[i-1];
01518 }
01519
01520 kints = fNiofex[ke1-1];
01521 mnfixp(kints-1, ierr);
01522 kints = fNiofex[ke2-1];
01523 mnfixp(kints-1, ierr);
01524
01525 for (inew = next; inew <= nptu; ++inew) {
01526
01527 bigdis = 0;
01528 for (iold = 1; iold <= inew - 1; ++iold) {
01529 i2 = iold + 1;
01530 if (i2 == inew) i2 = 1;
01531 d__1 = scalx*(xptu[iold-1] - xptu[i2-1]);
01532 d__2 = scaly*(yptu[iold-1] - yptu[i2-1]);
01533 dist = d__1*d__1 + d__2*d__2;
01534 if (dist > bigdis) {
01535 bigdis = dist;
01536 idist = iold;
01537 }
01538 }
01539 i1 = idist;
01540 i2 = i1 + 1;
01541 if (i2 == inew) i2 = 1;
01542
01543 a1 = .5;
01544 a2 = .5;
01545 L300:
01546 fXmidcr = a1*xptu[i1-1] + a2*xptu[i2-1];
01547 fYmidcr = a1*yptu[i1-1] + a2*yptu[i2-1];
01548 xdir = yptu[i2-1] - yptu[i1-1];
01549 ydir = xptu[i1-1] - xptu[i2-1];
01550 sclfac = TMath::Max(TMath::Abs(xdir*scalx),TMath::Abs(ydir*scaly));
01551 fXdircr = xdir / sclfac;
01552 fYdircr = ydir / sclfac;
01553 fKe1cr = ke1;
01554 fKe2cr = ke2;
01555
01556 fAmin = abest;
01557 mncros(aopt, iercr);
01558 if (iercr > 1) {
01559
01560 if (a1 > .5) {
01561 if (fISW[4] >= 0) {
01562 Printf(" MNCONT CANNOT FIND NEXT POINT ON CONTOUR. ONLY %3d POINTS FOUND.",nowpts);
01563 }
01564 goto L950;
01565 }
01566 mnwarn("W", "MNContour ", "Cannot find midpoint, try closer.");
01567 a1 = .75;
01568 a2 = .25;
01569 goto L300;
01570 }
01571
01572 for (move = nowpts; move >= i1 + 1; --move) {
01573 xptu[move] = xptu[move-1];
01574 yptu[move] = yptu[move-1];
01575 }
01576 ++nowpts;
01577 xptu[i1] = fXmidcr + fXdircr*aopt;
01578 yptu[i1] = fYmidcr + fYdircr*aopt;
01579 }
01580 L950:
01581
01582 ierrf = nowpts;
01583 fCstatu = "SUCCESSFUL";
01584 if (nowpts < nptu) fCstatu = "INCOMPLETE";
01585
01586
01587 if (fISW[4] >= 0) {
01588 fXpt[0] = u1min;
01589 fYpt[0] = u2min;
01590 fChpt[0] = ' ';
01591 nall = TMath::Min(nowpts + 1,101);
01592 for (i = 2; i <= nall; ++i) {
01593 fXpt[i-1] = xptu[i-2];
01594 fYpt[i-1] = yptu[i-2];
01595 fChpt[i-1] = 'X';
01596 }
01597 fChpt[nall] = 0;
01598 Printf(" Y-AXIS: PARAMETER %3d %s",ke2,(const char*)fCpnam[ke2-1]);
01599
01600 mnplot(fXpt, fYpt, fChpt, nall, fNpagwd, fNpagln);
01601
01602 Printf(" X-AXIS: PARAMETER %3d %s",ke1,(const char*)fCpnam[ke1-1]);
01603 }
01604
01605 if (fISW[4] >= 1) {
01606 npcol = (nowpts + 1) / 2;
01607 nfcol = nowpts / 2;
01608 Printf("%5d POINTS ON CONTOUR. FMIN=%13.5e ERRDEF=%11.3g",nowpts,abest,fUp);
01609 Printf(" %s%s%s%s",(const char*)fCpnam[ke1-1],
01610 (const char*)fCpnam[ke2-1],
01611 (const char*)fCpnam[ke1-1],
01612 (const char*)fCpnam[ke2-1]);
01613 for (line = 1; line <= nfcol; ++line) {
01614 lr = line + npcol;
01615 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]);
01616 }
01617 if (nfcol < npcol) {
01618 Printf(" %5d%13.5e%13.5e",npcol,xptu[npcol-1],yptu[npcol-1]);
01619 }
01620 }
01621
01622 fItaur = 1;
01623 mnfree(1);
01624 mnfree(1);
01625 i__1 = mpar*(mpar + 1) / 2;
01626 for (j = 1; j <= i__1; ++j) { fVhmat[j-1] = fVthmat[j-1]; }
01627 for (i = 1; i <= mpar; ++i) {
01628 fGlobcc[i-1] = fCONTgcc[i-1];
01629 fWerr[i-1] = fCONTw[i-1];
01630 fX[i-1] = fXt[i-1];
01631 }
01632 mninex(fX);
01633 fEDM = sigsav;
01634 fAmin = abest;
01635 fISW[1] = isw2;
01636 fISW[3] = isw4;
01637 fDcovar = dc;
01638 fItaur = 0;
01639 fNfcnmx = nfmxin;
01640 fIstrat = istrav;
01641 fU[ke1-1] = u1min;
01642 fU[ke2-1] = u2min;
01643 goto L2000;
01644
01645 L1350:
01646 Printf(" INVALID PARAMETER NUMBERS.");
01647 goto L1450;
01648 L1400:
01649 Printf(" LESS THAN FOUR POINTS REQUESTED.");
01650 L1450:
01651 ierrf = -1;
01652 fCstatu = "USER ERROR";
01653 goto L2000;
01654 L1500:
01655 Printf(" MNCONT UNABLE TO FIND FOUR POINTS.");
01656 fU[ke1-1] = u1min;
01657 fU[ke2-1] = u2min;
01658 ierrf = 0;
01659 fCstatu = "FAILED";
01660 L2000:
01661 fCfrom = "MNContour ";
01662 fNfcnfr = nfcnco;
01663 }
01664
01665
01666 void TMinuit::mncrck(TString cardbuf, Int_t maxcwd, TString &comand, Int_t &lnc,
01667 Int_t mxp, Double_t *plist, Int_t &llist, Int_t &ierr, Int_t)
01668 {
01669
01670
01671
01672
01673
01674
01675
01676
01677
01678
01679
01680
01681
01682 char *cnull = 0;
01683 const char *cnumer = "123456789-.0+";
01684
01685
01686 Int_t ifld, iend, lend, left, nreq, ipos, kcmnd, nextb, ic, ibegin, ltoadd;
01687 Int_t ielmnt, lelmnt[25], nelmnt;
01688 TString ctemp;
01689 char *celmnt[25];
01690 char command[25];
01691
01692
01693 char *crdbuf = (char*)cardbuf.Data();
01694 lend = cardbuf.Length();
01695 ielmnt = 0;
01696 nextb = 1;
01697 ierr = 0;
01698
01699 L10:
01700 for (ipos = nextb; ipos <= lend; ++ipos) {
01701 ibegin = ipos;
01702 if (crdbuf[ipos-1] == ' ') continue;
01703 if (crdbuf[ipos-1] == ',') goto L250;
01704 goto L150;
01705 }
01706 goto L300;
01707 L150:
01708
01709 for (ipos = ibegin + 1; ipos <= lend; ++ipos) {
01710 if (crdbuf[ipos-1] == ' ') goto L250;
01711 if (crdbuf[ipos-1] == ',') goto L250;
01712 }
01713 ipos = lend + 1;
01714 L250:
01715 iend = ipos - 1;
01716 ++ielmnt;
01717 if (iend >= ibegin) celmnt[ielmnt-1] = &crdbuf[ibegin-1];
01718 else celmnt[ielmnt-1] = cnull;
01719 lelmnt[ielmnt-1] = iend - ibegin + 1;
01720 if (lelmnt[ielmnt-1] > 19) {
01721 Printf(" MINUIT WARNING: INPUT DATA WORD TOO LONG.");
01722 ctemp = cardbuf(ibegin-1,iend-ibegin+1);
01723 Printf(" ORIGINAL:%s",ctemp.Data());
01724 Printf(" TRUNCATED TO:%s",celmnt[ielmnt-1]);
01725 lelmnt[ielmnt-1] = 19;
01726 }
01727 if (ipos >= lend) goto L300;
01728 if (ielmnt >= 25) goto L300;
01729
01730 for (ipos = iend + 1; ipos <= lend; ++ipos) {
01731 if (crdbuf[ipos-1] == ' ') continue;
01732 nextb = ipos;
01733 if (crdbuf[ipos-1] == ',') nextb = ipos + 1;
01734 goto L10;
01735 }
01736
01737
01738 L300:
01739 nelmnt = ielmnt;
01740 command[0] = ' '; command[1] = 0;
01741 lnc = 1;
01742 plist[0] = 0;
01743 llist = 0;
01744 if (ielmnt == 0) goto L900;
01745 kcmnd = 0;
01746 for (ielmnt = 1; ielmnt <= nelmnt; ++ielmnt) {
01747 if ( celmnt[ielmnt-1] == cnull) goto L450;
01748 for (ic = 1; ic <= 13; ++ic) {
01749 if (*celmnt[ielmnt-1] == cnumer[ic-1]) goto L450;
01750 }
01751 if (kcmnd >= maxcwd) continue;
01752 left = maxcwd - kcmnd;
01753 ltoadd = lelmnt[ielmnt-1];
01754 if (ltoadd > left) ltoadd = left;
01755 strncpy(&command[kcmnd],celmnt[ielmnt-1],ltoadd);
01756 kcmnd += ltoadd;
01757 if (kcmnd == maxcwd) continue;
01758 command[kcmnd] = ' ';
01759 ++kcmnd;
01760 command[kcmnd] = 0;
01761 }
01762 lnc = kcmnd;
01763 goto L900;
01764 L450:
01765 lnc = kcmnd;
01766
01767 llist = 0;
01768 for (ifld = ielmnt; ifld <= nelmnt; ++ifld) {
01769 ++llist;
01770 if (llist > mxp) {
01771 nreq = nelmnt - ielmnt + 1;
01772 Printf(" MINUIT WARNING IN MNCRCK: ");
01773 Printf(" COMMAND HAS INPUT %5d NUMERIC FIELDS, BUT MINUIT CAN ACCEPT ONLY%3d",nreq,mxp);
01774 goto L900;
01775 }
01776 if (celmnt[ifld-1] == cnull) plist[llist-1] = 0;
01777 else {
01778 sscanf(celmnt[ifld-1],"%lf",&plist[llist-1]);
01779 }
01780 }
01781
01782 L900:
01783 if (lnc <= 0) lnc = 1;
01784 comand = command;
01785 }
01786
01787
01788 void TMinuit::mncros(Double_t &aopt, Int_t &iercr)
01789 {
01790
01791
01792
01793
01794
01795
01796
01797
01798
01799
01800
01801
01802 Double_t alsb[3], flsb[3], bmin, bmax, zmid, sdev, zdir, zlim;
01803 Double_t coeff[3], aleft, aulim, fdist, adist, aminsv;
01804 Double_t anext, fnext, slope, s1, s2, x1, x2, ecarmn, ecarmx;
01805 Double_t determ, rt, smalla, aright, aim, tla, tlf, dfda,ecart;
01806 Int_t iout=0, i, ileft, ierev, maxlk, ibest, ik, it;
01807 Int_t noless, iworst=0, iright, itoohi, kex, ipt;
01808 Bool_t ldebug;
01809 const char *chsign;
01810 x2 = 0;
01811
01812 ldebug = fIdbg[6] >= 1;
01813 aminsv = fAmin;
01814
01815
01816 aim = fAmin + fUp;
01817 tlf = fUp*.01;
01818 tla = .01;
01819 fXpt[0] = 0;
01820 fYpt[0] = aim;
01821 fChpt[0] = ' ';
01822 ipt = 1;
01823 if (fKe2cr == 0) {
01824 fXpt[1] = -1;
01825 fYpt[1] = fAmin;
01826 fChpt[1] = '.';
01827 ipt = 2;
01828 }
01829
01830 aulim = 100;
01831 for (ik = 1; ik <= 2; ++ik) {
01832 if (ik == 1) {
01833 kex = fKe1cr;
01834 zmid = fXmidcr;
01835 zdir = fXdircr;
01836 } else {
01837 if (fKe2cr == 0) continue;
01838 kex = fKe2cr;
01839 zmid = fYmidcr;
01840 zdir = fYdircr;
01841 }
01842 if (fNvarl[kex-1] <= 1) continue;
01843 if (zdir == 0) continue;
01844 zlim = fAlim[kex-1];
01845 if (zdir > 0) zlim = fBlim[kex-1];
01846 aulim = TMath::Min(aulim,(zlim - zmid) / zdir);
01847 }
01848
01849
01850 anext = 0;
01851 aopt = anext;
01852 fLimset = kFALSE;
01853 if (aulim < aopt + tla) fLimset = kTRUE;
01854 mneval(anext, fnext, ierev);
01855
01856 if (ldebug) {
01857 Printf(" MNCROS: calls=%8d AIM=%10.5f F,A=%10.5f%10.5f",fNfcn,aim,fnext,aopt);
01858 }
01859 if (ierev > 0) goto L900;
01860 if (fLimset && fnext <= aim) goto L930;
01861 ++ipt;
01862 fXpt[ipt-1] = anext;
01863 fYpt[ipt-1] = fnext;
01864 fChpt[ipt-1] = charal[ipt-1];
01865 alsb[0] = anext;
01866 flsb[0] = fnext;
01867 fnext = TMath::Max(fnext,aminsv + fUp*.1);
01868 aopt = TMath::Sqrt(fUp / (fnext - aminsv)) - 1;
01869 if (TMath::Abs(fnext - aim) < tlf) goto L800;
01870
01871 if (aopt < -.5)aopt = -.5;
01872 if (aopt > 1) aopt = 1;
01873 fLimset = kFALSE;
01874 if (aopt > aulim) {
01875 aopt = aulim;
01876 fLimset = kTRUE;
01877 }
01878 mneval(aopt, fnext, ierev);
01879
01880 if (ldebug) {
01881 Printf(" MNCROS: calls=%8d AIM=%10.5f F,A=%10.5f%10.5f",fNfcn,aim,fnext,aopt);
01882 }
01883 if (ierev > 0) goto L900;
01884 if (fLimset && fnext <= aim) goto L930;
01885 alsb[1] = aopt;
01886 ++ipt;
01887 fXpt[ipt-1] = alsb[1];
01888 fYpt[ipt-1] = fnext;
01889 fChpt[ipt-1] = charal[ipt-1];
01890 flsb[1] = fnext;
01891 dfda = (flsb[1] - flsb[0]) / (alsb[1] - alsb[0]);
01892
01893 if (dfda > 0) goto L460;
01894 L300:
01895 mnwarn("D", "MNCROS ", "Looking for slope of the right sign");
01896 maxlk = 15 - ipt;
01897 for (it = 1; it <= maxlk; ++it) {
01898 alsb[0] = alsb[1];
01899 flsb[0] = flsb[1];
01900 aopt = alsb[0] + Double_t(it)*.2;
01901 fLimset = kFALSE;
01902 if (aopt > aulim) {
01903 aopt = aulim;
01904 fLimset = kTRUE;
01905 }
01906 mneval(aopt, fnext, ierev);
01907
01908 if (ldebug) {
01909 Printf(" MNCROS: calls=%8d AIM=%10.5f F,A=%10.5f%10.5f",fNfcn,aim,fnext,aopt);
01910 }
01911 if (ierev > 0) goto L900;
01912 if (fLimset && fnext <= aim) goto L930;
01913 alsb[1] = aopt;
01914 ++ipt;
01915 fXpt[ipt-1] = alsb[1];
01916 fYpt[ipt-1] = fnext;
01917 fChpt[ipt-1] = charal[ipt-1];
01918 flsb[1] = fnext;
01919 dfda = (flsb[1] - flsb[0]) / (alsb[1] - alsb[0]);
01920 if (dfda > 0) goto L450;
01921 }
01922 mnwarn("W", "MNCROS ", "Cannot find slope of the right sign");
01923 goto L950;
01924 L450:
01925
01926 L460:
01927 aopt = alsb[1] + (aim - flsb[1]) / dfda;
01928 fdist = TMath::Min(TMath::Abs(aim - flsb[0]),TMath::Abs(aim - flsb[1]));
01929 adist = TMath::Min(TMath::Abs(aopt - alsb[0]),TMath::Abs(aopt - alsb[1]));
01930 tla = .01;
01931 if (TMath::Abs(aopt) > 1) tla = TMath::Abs(aopt)*.01;
01932 if (adist < tla && fdist < tlf) goto L800;
01933 if (ipt >= 15) goto L950;
01934 bmin = TMath::Min(alsb[0],alsb[1]) - 1;
01935 if (aopt < bmin) aopt = bmin;
01936 bmax = TMath::Max(alsb[0],alsb[1]) + 1;
01937 if (aopt > bmax) aopt = bmax;
01938
01939 fLimset = kFALSE;
01940 if (aopt > aulim) {
01941 aopt = aulim;
01942 fLimset = kTRUE;
01943 }
01944 mneval(aopt, fnext, ierev);
01945
01946 if (ldebug) {
01947 Printf(" MNCROS: calls=%8d AIM=%10.5f F,A=%10.5f%10.5f",fNfcn,aim,fnext,aopt);
01948 }
01949 if (ierev > 0) goto L900;
01950 if (fLimset && fnext <= aim) goto L930;
01951 alsb[2] = aopt;
01952 ++ipt;
01953 fXpt[ipt-1] = alsb[2];
01954 fYpt[ipt-1] = fnext;
01955 fChpt[ipt-1] = charal[ipt-1];
01956 flsb[2] = fnext;
01957
01958 ecarmn = TMath::Abs(fnext-aim);
01959 ibest = 3;
01960 ecarmx = 0;
01961 noless = 0;
01962 for (i = 1; i <= 3; ++i) {
01963 ecart = TMath::Abs(flsb[i-1] - aim);
01964 if (ecart > ecarmx) { ecarmx = ecart; iworst = i; }
01965 if (ecart < ecarmn) { ecarmn = ecart; ibest = i; }
01966 if (flsb[i-1] < aim) ++noless;
01967 }
01968
01969 if (noless == 1 || noless == 2) goto L500;
01970
01971 if (noless == 0 && ibest != 3) goto L950;
01972
01973
01974 if (noless == 3 && ibest != 3) {
01975 alsb[1] = alsb[2];
01976 flsb[1] = flsb[2];
01977 goto L300;
01978 }
01979
01980 alsb[iworst-1] = alsb[2];
01981 flsb[iworst-1] = flsb[2];
01982 dfda = (flsb[1] - flsb[0]) / (alsb[1] - alsb[0]);
01983 goto L460;
01984
01985 L500:
01986 mnpfit(alsb, flsb, 3, coeff, sdev);
01987 if (coeff[2] <= 0) {
01988 mnwarn("D", "MNCROS ", "Curvature is negative near contour line.");
01989 }
01990 determ = coeff[1]*coeff[1] - coeff[2]*4*(coeff[0] - aim);
01991 if (determ <= 0) {
01992 mnwarn("D", "MNCROS ", "Problem 2, impossible determinant");
01993 goto L950;
01994 }
01995
01996 rt = TMath::Sqrt(determ);
01997 x1 = (-coeff[1] + rt) / (coeff[2]*2);
01998 x2 = (-coeff[1] - rt) / (coeff[2]*2);
01999 s1 = coeff[1] + x1*2*coeff[2];
02000 s2 = coeff[1] + x2*2*coeff[2];
02001 if (s1*s2 > 0) {
02002 Printf(" MNCONTour problem 1");
02003 }
02004 aopt = x1;
02005 slope = s1;
02006 if (s2 > 0) {
02007 aopt = x2;
02008 slope = s2;
02009 }
02010
02011 tla = .01;
02012 if (TMath::Abs(aopt) > 1) tla = TMath::Abs(aopt)*.01;
02013 if (TMath::Abs(aopt - alsb[ibest-1]) < tla && TMath::Abs(flsb[ibest-1] - aim) < tlf) {
02014 goto L800;
02015 }
02016 if (ipt >= 15) goto L950;
02017
02018
02019
02020 ileft = 0;
02021 iright = 0;
02022 ibest = 1;
02023 ecarmx = 0;
02024 ecarmn = TMath::Abs(aim - flsb[0]);
02025 for (i = 1; i <= 3; ++i) {
02026 ecart = TMath::Abs(flsb[i-1] - aim);
02027 if (ecart < ecarmn) { ecarmn = ecart; ibest = i; }
02028 if (ecart > ecarmx) { ecarmx = ecart; }
02029 if (flsb[i-1] > aim) {
02030 if (iright == 0) iright = i;
02031 else if (flsb[i-1] > flsb[iright-1]) iout = i;
02032 else { iout = iright; iright = i; }
02033 }
02034 else if (ileft == 0) ileft = i;
02035 else if (flsb[i-1] < flsb[ileft-1]) iout = i;
02036 else { iout = ileft; ileft = i; }
02037 }
02038
02039 if (ecarmx > TMath::Abs(flsb[iout-1] - aim)*10) {
02040 aopt = aopt*.5 + (alsb[iright-1] + alsb[ileft-1])*.25;
02041 }
02042
02043 smalla = tla*.1;
02044 if (slope*smalla > tlf) smalla = tlf / slope;
02045 aleft = alsb[ileft-1] + smalla;
02046 aright = alsb[iright-1] - smalla;
02047
02048 if (aopt < aleft) aopt = aleft;
02049 if (aopt > aright) aopt = aright;
02050 if (aleft > aright) aopt = (aleft + aright)*.5;
02051
02052
02053 fLimset = kFALSE;
02054 if (aopt > aulim) {
02055 aopt = aulim;
02056 fLimset = kTRUE;
02057 }
02058
02059 mneval(aopt, fnext, ierev);
02060
02061 if (ldebug) {
02062 Printf(" MNCROS: calls=%8d AIM=%10.5f F,A=%10.5f%10.5f",fNfcn,aim,fnext,aopt);
02063 }
02064 if (ierev > 0) goto L900;
02065 if (fLimset && fnext <= aim) goto L930;
02066 ++ipt;
02067 fXpt[ipt-1] = aopt;
02068 fYpt[ipt-1] = fnext;
02069 fChpt[ipt-1] = charal[ipt-1];
02070
02071 alsb[iout-1] = aopt;
02072 flsb[iout-1] = fnext;
02073
02074
02075 ibest = iout;
02076 goto L500;
02077
02078
02079 L800:
02080 iercr = 0;
02081 goto L1000;
02082
02083 L900:
02084 if (ierev == 1) goto L940;
02085 goto L950;
02086
02087 L930:
02088 iercr = 1;
02089 goto L1000;
02090
02091 L940:
02092 iercr = 2;
02093 goto L1000;
02094
02095 L950:
02096 iercr = 3;
02097
02098 L1000:
02099 if (ldebug) {
02100 itoohi = 0;
02101 for (i = 1; i <= ipt; ++i) {
02102 if (fYpt[i-1] > aim + fUp) {
02103 fYpt[i-1] = aim + fUp;
02104 fChpt[i-1] = '+';
02105 itoohi = 1;
02106 }
02107 }
02108 fChpt[ipt] = 0;
02109 chsign = "POSI";
02110 if (fXdircr < 0) chsign = "NEGA";
02111 if (fKe2cr == 0) {
02112 Printf(" %sTIVE MINOS ERROR, PARAMETER %3d",chsign,fKe1cr);
02113 }
02114 if (itoohi == 1) {
02115 Printf("POINTS LABELLED '+' WERE TOO HIGH TO PLOT.");
02116 }
02117 if (iercr == 1) {
02118 Printf("RIGHTMOST POINT IS UP AGAINST LIMIT.");
02119 }
02120 mnplot(fXpt, fYpt, fChpt, ipt, fNpagwd, fNpagln);
02121 }
02122 }
02123
02124
02125 void TMinuit::mncuve()
02126 {
02127
02128
02129
02130
02131
02132
02133
02134
02135 Double_t dxdi, wint;
02136 Int_t ndex, iext, i, j;
02137
02138 if (fISW[3] < 1) {
02139 Printf(" FUNCTION MUST BE MINIMIZED BEFORE CALLING %s",(const char*)fCfrom);
02140 fApsi = fEpsi;
02141 mnmigr();
02142 }
02143 if (fISW[1] < 3) {
02144 mnhess();
02145 if (fISW[1] < 1) {
02146 mnwarn("W", fCfrom, "NO ERROR MATRIX. WILL IMPROVISE.");
02147 for (i = 1; i <= fNpar; ++i) {
02148 ndex = i*(i-1) / 2;
02149 for (j = 1; j <= i-1; ++j) {
02150 ++ndex;
02151 fVhmat[ndex-1] = 0;
02152 }
02153 ++ndex;
02154 if (fG2[i-1] <= 0) {
02155 wint = fWerr[i-1];
02156 iext = fNexofi[i-1];
02157 if (fNvarl[iext-1] > 1) {
02158 mndxdi(fX[i-1], i-1, dxdi);
02159 if (TMath::Abs(dxdi) < .001) wint = .01;
02160 else wint /= TMath::Abs(dxdi);
02161 }
02162 fG2[i-1] = fUp / (wint*wint);
02163 }
02164 fVhmat[ndex-1] = 2 / fG2[i-1];
02165 }
02166 fISW[1] = 1;
02167 fDcovar = 1;
02168 } else mnwerr();
02169 }
02170 }
02171
02172
02173 void TMinuit::mnderi()
02174 {
02175
02176
02177
02178
02179
02180
02181
02182
02183
02184 Double_t step, dfmin, stepb4, dd, df, fs1;
02185 Double_t tlrstp, tlrgrd, epspri, optstp, stpmax, stpmin, fs2, grbfor=0, d1d2, xtf;
02186 Int_t icyc, ncyc, iint, iext, i, nparx;
02187 Bool_t ldebug;
02188
02189 nparx = fNpar;
02190 ldebug = fIdbg[2] >= 1;
02191 if (fAmin == fUndefi) mnamin();
02192 if (fISW[2] == 1) goto L100;
02193
02194 if (ldebug) {
02195
02196 mninex(fX);
02197 nparx = fNpar;
02198 Eval(nparx, fGin, fs1, fU, 4); ++fNfcn;
02199 if (fs1 != fAmin) {
02200 df = fAmin - fs1;
02201 mnwarn("D", "MNDERI", Form("function value differs from AMIN by %12.3g",df));
02202 fAmin = fs1;
02203 }
02204 Printf(" FIRST DERIVATIVE DEBUG PRINTOUT. MNDERI");
02205 Printf(" PAR DERIV STEP MINSTEP OPTSTEP D1-D2 2ND DRV");
02206 }
02207 dfmin = fEpsma2*8*(TMath::Abs(fAmin) + fUp);
02208 if (fIstrat <= 0) {
02209 ncyc = 2;
02210 tlrstp = .5;
02211 tlrgrd = .1;
02212 } else if (fIstrat == 1) {
02213 ncyc = 3;
02214 tlrstp = .3;
02215 tlrgrd = .05;
02216 } else {
02217 ncyc = 5;
02218 tlrstp = .1;
02219 tlrgrd = .02;
02220 }
02221
02222 for (i = 1; i <= fNpar; ++i) {
02223 epspri = fEpsma2 + TMath::Abs(fGrd[i-1]*fEpsma2);
02224
02225
02226 xtf = fX[i-1];
02227 stepb4 = 0;
02228
02229 for (icyc = 1; icyc <= ncyc; ++icyc) {
02230
02231 optstp = TMath::Sqrt(dfmin / (TMath::Abs(fG2[i-1]) + epspri));
02232
02233 step = TMath::Max(optstp,TMath::Abs(fGstep[i-1]*.1));
02234
02235 if (fGstep[i-1] < 0 && step > .5) step = .5;
02236
02237 stpmax = TMath::Abs(fGstep[i-1])*10;
02238 if (step > stpmax) step = stpmax;
02239
02240 stpmin = TMath::Abs(fEpsma2*fX[i-1])*8;
02241 if (step < stpmin) step = stpmin;
02242
02243 if (TMath::Abs((step - stepb4) / step) < tlrstp) goto L50;
02244
02245 stepb4 = step;
02246 if (fGstep[i-1] > 0) fGstep[i-1] = TMath::Abs(step);
02247 else fGstep[i-1] = -TMath::Abs(step);
02248 stepb4 = step;
02249 fX[i-1] = xtf + step;
02250 mninex(fX);
02251 Eval(nparx, fGin, fs1, fU, 4); ++fNfcn;
02252
02253 fX[i-1] = xtf - step;
02254 mninex(fX);
02255 Eval(nparx, fGin, fs2, fU, 4); ++fNfcn;
02256 grbfor = fGrd[i-1];
02257 fGrd[i-1] = (fs1 - fs2) / (step*2);
02258 fG2[i-1] = (fs1 + fs2 - fAmin*2) / (step*step);
02259 fX[i-1] = xtf;
02260 if (ldebug) {
02261 d1d2 = (fs1 + fs2 - fAmin*2) / step;
02262 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]);
02263 }
02264
02265 if (TMath::Abs(grbfor - fGrd[i-1]) / (TMath::Abs(fGrd[i-1]) + dfmin/step) < tlrgrd)
02266 goto L50;
02267 }
02268
02269 if (ncyc == 1) goto L50;
02270 mnwarn("D", "MNDERI", Form("First derivative not converged. %g%g",fGrd[i-1],grbfor));
02271 L50:
02272 ;
02273 }
02274 mninex(fX);
02275 return;
02276
02277 L100:
02278 for (iint = 1; iint <= fNpar; ++iint) {
02279 iext = fNexofi[iint-1];
02280 if (fNvarl[iext-1] <= 1) {
02281 fGrd[iint-1] = fGin[iext-1];
02282 } else {
02283 dd = (fBlim[iext-1] - fAlim[iext-1])*.5*TMath::Cos(fX[iint-1]);
02284 fGrd[iint-1] = fGin[iext-1]*dd;
02285 }
02286 }
02287 }
02288
02289
02290 void TMinuit::mndxdi(Double_t pint, Int_t ipar, Double_t &dxdi)
02291 {
02292
02293
02294
02295
02296
02297
02298
02299 Int_t i = fNexofi[ipar];
02300 dxdi = 1;
02301 if (fNvarl[i-1] > 1) {
02302 dxdi = TMath::Abs((fBlim[i-1] - fAlim[i-1])*TMath::Cos(pint))*.5;
02303 }
02304 }
02305
02306
02307 void TMinuit::mneig(Double_t *a, Int_t ndima, Int_t n, Int_t mits, Double_t *work, Double_t precis, Int_t &ifault)
02308 {
02309
02310
02311
02312 Int_t a_offset;
02313 Double_t d__1;
02314
02315
02316 Double_t b, c, f, h, r, s, hh, gl, pr, pt;
02317 Int_t i, j, k, l, m=0, i0, i1, j1, m1, n1;
02318
02319
02320
02321 a_offset = ndima + 1;
02322 a -= a_offset;
02323 --work;
02324
02325
02326 ifault = 1;
02327
02328 i = n;
02329 for (i1 = 2; i1 <= n; ++i1) {
02330 l = i-2;
02331 f = a[i + (i-1)*ndima];
02332 gl = 0;
02333
02334 if (l < 1) goto L25;
02335
02336 for (k = 1; k <= l; ++k) {
02337 d__1 = a[i + k*ndima];
02338 gl += d__1*d__1;
02339 }
02340 L25:
02341 h = gl + f*f;
02342
02343 if (gl > 1e-35) goto L30;
02344
02345 work[i] = 0;
02346 work[n + i] = f;
02347 goto L65;
02348 L30:
02349 ++l;
02350 gl = TMath::Sqrt(h);
02351 if (f >= 0) gl = -gl;
02352 work[n + i] = gl;
02353 h -= f*gl;
02354 a[i + (i-1)*ndima] = f - gl;
02355 f = 0;
02356 for (j = 1; j <= l; ++j) {
02357 a[j + i*ndima] = a[i + j*ndima] / h;
02358 gl = 0;
02359 for (k = 1; k <= j; ++k) { gl += a[j + k*ndima]*a[i + k*ndima]; }
02360 if (j >= l) goto L47;
02361 j1 = j + 1;
02362 for (k = j1; k <= l; ++k) { gl += a[k + j*ndima]*a[i + k*ndima]; }
02363 L47:
02364 work[n + j] = gl / h;
02365 f += gl*a[j + i*ndima];
02366 }
02367 hh = f / (h + h);
02368 for (j = 1; j <= l; ++j) {
02369 f = a[i + j*ndima];
02370 gl = work[n + j] - hh*f;
02371 work[n + j] = gl;
02372 for (k = 1; k <= j; ++k) {
02373 a[j + k*ndima] = a[j + k*ndima] - f*work[n + k] - gl*a[i + k*ndima];
02374 }
02375 }
02376 work[i] = h;
02377 L65:
02378 --i;
02379 }
02380 work[1] = 0;
02381 work[n + 1] = 0;
02382 for (i = 1; i <= n; ++i) {
02383 l = i-1;
02384 if (work[i] == 0 || l == 0) goto L100;
02385
02386 for (j = 1; j <= l; ++j) {
02387 gl = 0;
02388 for (k = 1; k <= l; ++k) { gl += a[i + k*ndima]*a[k + j*ndima]; }
02389 for (k = 1; k <= l; ++k) { a[k + j*ndima] -= gl*a[k + i*ndima]; }
02390 }
02391 L100:
02392 work[i] = a[i + i*ndima];
02393 a[i + i*ndima] = 1;
02394 if (l == 0) continue;
02395
02396 for (j = 1; j <= l; ++j) {
02397 a[i + j*ndima] = 0;
02398 a[j + i*ndima] = 0;
02399 }
02400 }
02401
02402 n1 = n - 1;
02403 for (i = 2; i <= n; ++i) {
02404 i0 = n + i-1;
02405 work[i0] = work[i0 + 1];
02406 }
02407 work[n + n] = 0;
02408 b = 0;
02409 f = 0;
02410 for (l = 1; l <= n; ++l) {
02411 j = 0;
02412 h = precis*(TMath::Abs(work[l]) + TMath::Abs(work[n + l]));
02413 if (b < h) b = h;
02414 for (m1 = l; m1 <= n; ++m1) {
02415 m = m1;
02416 if (TMath::Abs(work[n + m]) <= b) goto L150;
02417 }
02418
02419 L150:
02420 if (m == l) goto L205;
02421
02422 L160:
02423 if (j == mits) return;
02424 ++j;
02425 pt = (work[l + 1] - work[l]) / (work[n + l]*2);
02426 r = TMath::Sqrt(pt*pt + 1);
02427 pr = pt + r;
02428 if (pt < 0) pr = pt - r;
02429
02430 h = work[l] - work[n + l] / pr;
02431 for (i = l; i <= n; ++i) { work[i] -= h; }
02432 f += h;
02433 pt = work[m];
02434 c = 1;
02435 s = 0;
02436 m1 = m - 1;
02437 i = m;
02438 for (i1 = l; i1 <= m1; ++i1) {
02439 j = i;
02440 --i;
02441 gl = c*work[n + i];
02442 h = c*pt;
02443 if (TMath::Abs(pt) >= TMath::Abs(work[n + i])) goto L180;
02444
02445 c = pt / work[n + i];
02446 r = TMath::Sqrt(c*c + 1);
02447 work[n + j] = s*work[n + i]*r;
02448 s = 1 / r;
02449 c /= r;
02450 goto L190;
02451 L180:
02452 c = work[n + i] / pt;
02453 r = TMath::Sqrt(c*c + 1);
02454 work[n + j] = s*pt*r;
02455 s = c / r;
02456 c = 1 / r;
02457 L190:
02458 pt = c*work[i] - s*gl;
02459 work[j] = h + s*(c*gl + s*work[i]);
02460 for (k = 1; k <= n; ++k) {
02461 h = a[k + j*ndima];
02462 a[k + j*ndima] = s*a[k + i*ndima] + c*h;
02463 a[k + i*ndima] = c*a[k + i*ndima] - s*h;
02464 }
02465 }
02466 work[n + l] = s*pt;
02467 work[l] = c*pt;
02468
02469 if (TMath::Abs(work[n + l]) > b) goto L160;
02470
02471 L205:
02472 work[l] += f;
02473 }
02474 for (i = 1; i <= n1; ++i) {
02475 k = i;
02476 pt = work[i];
02477 i1 = i + 1;
02478 for (j = i1; j <= n; ++j) {
02479 if (work[j] >= pt) continue;
02480 k = j;
02481 pt = work[j];
02482 }
02483
02484 if (k == i) continue;
02485
02486 work[k] = work[i];
02487 work[i] = pt;
02488 for (j = 1; j <= n; ++j) {
02489 pt = a[j + i*ndima];
02490 a[j + i*ndima] = a[j + k*ndima];
02491 a[j + k*ndima] = pt;
02492 }
02493 }
02494 ifault = 0;
02495 }
02496
02497
02498 void TMinuit::mnemat(Double_t *emat, Int_t ndim)
02499 {
02500
02501
02502
02503
02504
02505
02506
02507 Int_t emat_dim1, emat_offset;
02508
02509
02510 Double_t dxdi, dxdj;
02511 Int_t i, j, k, npard, k2, kk, iz, nperln, kga, kgb;
02512 TString ctemp;
02513
02514
02515 emat_dim1 = ndim;
02516 emat_offset = emat_dim1 + 1;
02517 emat -= emat_offset;
02518
02519
02520 if (fISW[1] < 1) return;
02521 if (fISW[4] >= 2) {
02522 Printf(" EXTERNAL ERROR MATRIX. NDIM=%4d NPAR=%3d ERR DEF=%g",ndim,fNpar,fUp);
02523 }
02524
02525 npard = fNpar;
02526 if (ndim < fNpar) {
02527 npard = ndim;
02528 if (fISW[4] >= 0) {
02529 Printf(" USER-DIMENSIONED ARRAY EMAT NOT BIG ENOUGH. REDUCED MATRIX CALCULATED.");
02530 }
02531 }
02532
02533
02534 nperln = (fNpagwd - 5) / 10;
02535 nperln = TMath::Min(nperln,13);
02536 if (fISW[4] >= 1 && npard > nperln) {
02537 Printf(" ELEMENTS ABOVE DIAGONAL ARE NOT PRINTED.");
02538 }
02539
02540 for (i = 1; i <= npard; ++i) {
02541 mndxdi(fX[i-1], i-1, dxdi);
02542 kga = i*(i-1) / 2;
02543 for (j = 1; j <= i; ++j) {
02544 mndxdi(fX[j-1], j-1, dxdj);
02545 kgb = kga + j;
02546 emat[i + j*emat_dim1] = dxdi*fVhmat[kgb-1]*dxdj*fUp;
02547 emat[j + i*emat_dim1] = emat[i + j*emat_dim1];
02548 }
02549 }
02550
02551 if (fISW[4] >= 2) {
02552 for (i = 1; i <= npard; ++i) {
02553 iz = npard;
02554 if (npard >= nperln) iz = i;
02555 ctemp = " ";
02556 for (k = 1; nperln < 0 ? k >= iz : k <= iz; k += nperln) {
02557 k2 = k + nperln - 1;
02558 if (k2 > iz) k2 = iz;
02559 for (kk = k; kk <= k2; ++kk) {
02560 ctemp += Form("%10.3e ",emat[i + kk*emat_dim1]);
02561 }
02562 Printf("%s",(const char*)ctemp);
02563 }
02564 }
02565 }
02566 }
02567
02568
02569 void TMinuit::mnerrs(Int_t number, Double_t &eplus, Double_t &eminus, Double_t &eparab, Double_t &gcc)
02570 {
02571
02572
02573
02574
02575
02576
02577
02578
02579
02580
02581
02582 Double_t dxdi;
02583 Int_t ndiag, iin, iex;
02584
02585 iex = number+1;
02586
02587 if (iex > fNu || iex <= 0) goto L900;
02588 iin = fNiofex[iex-1];
02589 if (iin <= 0) goto L900;
02590
02591
02592 eplus = fErp[iin-1];
02593 if (eplus == fUndefi) eplus = 0;
02594 eminus = fErn[iin-1];
02595 if (eminus == fUndefi) eminus = 0;
02596 mndxdi(fX[iin-1], iin-1, dxdi);
02597 ndiag = iin*(iin + 1) / 2;
02598 eparab = TMath::Abs(dxdi*TMath::Sqrt(TMath::Abs(fUp*fVhmat[ndiag- 1])));
02599
02600 gcc = 0;
02601 if (fISW[1] < 2) return;
02602 gcc = fGlobcc[iin-1];
02603 return;
02604
02605 L900:
02606 eplus = 0;
02607 eminus = 0;
02608 eparab = 0;
02609 gcc = 0;
02610 }
02611
02612
02613 void TMinuit::mneval(Double_t anext, Double_t &fnext, Int_t &ierev)
02614 {
02615
02616
02617
02618
02619
02620
02621
02622
02623
02624
02625 Int_t nparx;
02626
02627 fU[fKe1cr-1] = fXmidcr + anext*fXdircr;
02628 if (fKe2cr != 0) fU[fKe2cr-1] = fYmidcr + anext*fYdircr;
02629 mninex(fX);
02630 nparx = fNpar;
02631 Eval(nparx, fGin, fnext, fU, 4); ++fNfcn;
02632 ierev = 0;
02633 if (fNpar > 0) {
02634 fItaur = 1;
02635 fAmin = fnext;
02636 fISW[0] = 0;
02637 mnmigr();
02638 fItaur = 0;
02639 fnext = fAmin;
02640 if (fISW[0] >= 1) ierev = 1;
02641 if (fISW[3] < 1) ierev = 2;
02642 }
02643 }
02644
02645
02646 void TMinuit::mnexcm(const char *command, Double_t *plist, Int_t llist, Int_t &ierflg)
02647 {
02648
02649
02650
02651
02652
02653
02654
02655
02656
02657
02658
02659
02660
02661
02662
02663
02664
02665
02666
02667
02668
02669
02670
02671
02672
02673 TString comand = command;
02674 static TString clower = "abcdefghijklmnopqrstuvwxyz";
02675 static TString cupper = "ABCDEFGHIJKLMNOPQRSTUVWXYZ";
02676 const char *cname[40] = {
02677 "MINImize ",
02678 "SEEk ",
02679 "SIMplex ",
02680 "MIGrad ",
02681 "MINOs ",
02682 "SET xxx ",
02683 "SHOw xxx ",
02684 "TOP of pag",
02685 "FIX ",
02686 "REStore ",
02687 "RELease ",
02688 "SCAn ",
02689 "CONtour ",
02690 "HESse ",
02691 "SAVe ",
02692 "IMProve ",
02693 "CALl fcn ",
02694 "STAndard ",
02695 "END ",
02696 "EXIt ",
02697 "RETurn ",
02698 "CLEar ",
02699 "HELP ",
02700 "MNContour ",
02701 "STOp ",
02702 "JUMp ",
02703 " ",
02704 " ",
02705 " ",
02706 " ",
02707 " ",
02708 " ",
02709 " ",
02710 "COVARIANCE",
02711 "PRINTOUT ",
02712 "GRADIENT ",
02713 "MATOUT ",
02714 "ERROR DEF ",
02715 "LIMITS ",
02716 "PUNCH "};
02717
02718 Int_t nntot = 40;
02719
02720
02721 Double_t step, xptu[101], yptu[101], f, rno;
02722 Int_t icol, kcol, ierr, iint, iext, lnow, nptu, i, iflag, ierrf;
02723 Int_t ilist, nparx, izero, nf, lk, it, iw, inonde, nsuper;
02724 Int_t it2, ke1, ke2, nowprt, kll, krl;
02725 TString chwhy, c26, cvblnk, cneway, comd;
02726 TString ctemp;
02727 Bool_t lfreed, ltofix, lfixed;
02728
02729
02730
02731
02732
02733 lk = comand.Length();
02734 if (lk > 20) lk = 20;
02735 fCword = comand;
02736 fCword.ToUpper();
02737
02738
02739 for (iw = 1; iw <= fMaxpar; ++iw) {
02740 fWord7[iw-1] = 0;
02741 if (iw <= llist) fWord7[iw-1] = plist[iw-1];
02742 }
02743 ++fIcomnd;
02744 fNfcnlc = fNfcn;
02745 if (fCword(0,7) != "SET PRI" || fWord7[0] >= 0) {
02746 if (fISW[4] >= 0) {
02747 lnow = llist;
02748 if (lnow > 4) lnow = 4;
02749 Printf(" **********");
02750 ctemp = Form(" **%5d **%s",fIcomnd,(const char*)fCword);
02751 for (i = 1; i <= lnow; ++i) {
02752 ctemp += Form("%12.4g",plist[i-1]);
02753 }
02754 Printf("%s",(const char*)ctemp);
02755 inonde = 0;
02756 if (llist > lnow) {
02757 kll = llist;
02758 if (llist > fMaxpar) {
02759 inonde = 1;
02760 kll = fMaxpar;
02761 }
02762 Printf(" ***********");
02763 for (i = lnow + 1; i <= kll; ++i) {
02764 Printf("%12.4g",plist[i-1]);
02765 }
02766 }
02767 Printf(" **********");
02768 if (inonde > 0) {
02769 Printf(" ERROR: ABOVE CALL TO MNEXCM TRIED TO PASS MORE THAN %d PARAMETERS.", fMaxpar);
02770 }
02771 }
02772 }
02773 fNfcnmx = Int_t(fWord7[0]);
02774 if (fNfcnmx <= 0) {
02775 fNfcnmx = fNpar*100 + 200 + fNpar*fNpar*5;
02776 }
02777 fEpsi = fWord7[1];
02778 if (fEpsi <= 0) {
02779 fEpsi = fUp*.1;
02780 }
02781 fLnewmn = kFALSE;
02782 fLphead = kTRUE;
02783 fISW[0] = 0;
02784 ierflg = 0;
02785
02786 ctemp = fCword(0,3);
02787 for (i = 1; i <= nntot; ++i) {
02788 if (strncmp(ctemp.Data(),cname[i-1],3) == 0) goto L90;
02789 }
02790 Printf("UNKNOWN COMMAND IGNORED:%s", comand.Data());
02791 ierflg = 3;
02792 return;
02793
02794 L90:
02795 if (fCword(0,4) == "MINO") i = 5;
02796 if (i != 6 && i != 7 && i != 8 && i != 23) {
02797 fCfrom = cname[i-1];
02798 fNfcnfr = fNfcn;
02799 }
02800
02801 switch (i) {
02802 case 1: goto L400;
02803 case 2: goto L200;
02804 case 3: goto L300;
02805 case 4: goto L400;
02806 case 5: goto L500;
02807 case 6: goto L700;
02808 case 7: goto L700;
02809 case 8: goto L800;
02810 case 9: goto L900;
02811 case 10: goto L1000;
02812 case 11: goto L1100;
02813 case 12: goto L1200;
02814 case 13: goto L1300;
02815 case 14: goto L1400;
02816 case 15: goto L1500;
02817 case 16: goto L1600;
02818 case 17: goto L1700;
02819 case 18: goto L1800;
02820 case 19: goto L1900;
02821 case 20: goto L1900;
02822 case 21: goto L1900;
02823 case 22: goto L2200;
02824 case 23: goto L2300;
02825 case 24: goto L2400;
02826 case 25: goto L1900;
02827 case 26: goto L2600;
02828 case 27: goto L3300;
02829 case 28: goto L3300;
02830 case 29: goto L3300;
02831 case 30: goto L3300;
02832 case 31: goto L3300;
02833 case 32: goto L3300;
02834 case 33: goto L3300;
02835 case 34: goto L3400;
02836 case 35: goto L3500;
02837 case 36: goto L3600;
02838 case 37: goto L3700;
02839 case 38: goto L3800;
02840 case 39: goto L3900;
02841 case 40: goto L4000;
02842 }
02843
02844 L200:
02845 mnseek();
02846 return;
02847
02848 L300:
02849 mnsimp();
02850 if (fISW[3] < 1) ierflg = 4;
02851 return;
02852
02853 L400:
02854 nf = fNfcn;
02855 fApsi = fEpsi;
02856 mnmigr();
02857 mnwerr();
02858 if (fISW[3] >= 1) return;
02859 ierflg = 4;
02860 if (fISW[0] == 1) return;
02861 if (fCword(0,3) == "MIG") return;
02862
02863 fNfcnmx = fNfcnmx + nf - fNfcn;
02864 nf = fNfcn;
02865 mnsimp();
02866 if (fISW[0] == 1) return;
02867 fNfcnmx = fNfcnmx + nf - fNfcn;
02868 mnmigr();
02869 if (fISW[3] >= 1) ierflg = 0;
02870 mnwerr();
02871 return;
02872
02873 L500:
02874 nsuper = fNfcn + ((fNpar + 1) << 1)*fNfcnmx;
02875
02876 fEpsi = fUp*.1;
02877 L510:
02878 mncuve();
02879 mnmnos();
02880 if (! fLnewmn) return;
02881 mnrset(0);
02882 mnmigr();
02883 mnwerr();
02884 if (fNfcn < nsuper) goto L510;
02885 Printf(" TOO MANY FUNCTION CALLS. MINOS GIVES UP");
02886 ierflg = 4;
02887 return;
02888
02889 L700:
02890 mnset();
02891 return;
02892
02893
02894 L800:
02895 Printf("1");
02896 return;
02897
02898 L900:
02899 ltofix = kTRUE;
02900
02901 L901:
02902 lfreed = kFALSE;
02903 lfixed = kFALSE;
02904 if (llist == 0) {
02905 Printf("%s: NO PARAMETERS REQUESTED ",(const char*)fCword);
02906 return;
02907 }
02908 for (ilist = 1; ilist <= llist; ++ilist) {
02909 iext = Int_t(plist[ilist-1]);
02910 chwhy = " IS UNDEFINED.";
02911 if (iext <= 0) goto L930;
02912 if (iext > fNu) goto L930;
02913 if (fNvarl[iext-1] < 0) goto L930;
02914 chwhy = " IS CONSTANT. ";
02915 if (fNvarl[iext-1] == 0) goto L930;
02916 iint = fNiofex[iext-1];
02917 if (ltofix) {
02918 chwhy = " ALREADY FIXED.";
02919 if (iint == 0) goto L930;
02920 mnfixp(iint-1, ierr);
02921 if (ierr == 0) lfixed = kTRUE;
02922 else ierflg = 4;
02923 } else {
02924 chwhy = " ALREADY VARIABLE.";
02925 if (iint > 0) goto L930;
02926 krl = -abs(iext);
02927 mnfree(krl);
02928 lfreed = kTRUE;
02929 }
02930 continue;
02931 L930:
02932 if (fISW[4] >= 0) Printf(" PARAMETER %4d %s IGNORED.",iext,(const char*)chwhy);
02933 }
02934 if (lfreed || lfixed) mnrset(0);
02935 if (lfreed) {
02936 fISW[1] = 0;
02937 fDcovar = 1;
02938 fEDM = fBigedm;
02939 fISW[3] = 0;
02940 }
02941 mnwerr();
02942 if (fISW[4] > 1) mnprin(5, fAmin);
02943 return;
02944
02945 L1000:
02946 it = Int_t(fWord7[0]);
02947 if (it > 1 || it < 0) goto L1005;
02948 lfreed = fNpfix > 0;
02949 mnfree(it);
02950 if (lfreed) {
02951 mnrset(0);
02952 fISW[1] = 0;
02953 fDcovar = 1;
02954 fEDM = fBigedm;
02955 }
02956 return;
02957 L1005:
02958 Printf(" IGNORED. UNKNOWN ARGUMENT:%4d",it);
02959 ierflg = 3;
02960 return;
02961
02962 L1100:
02963 ltofix = kFALSE;
02964 goto L901;
02965
02966 L1200:
02967 iext = Int_t(fWord7[0]);
02968 if (iext <= 0) goto L1210;
02969 it2 = 0;
02970 if (iext <= fNu) it2 = fNiofex[iext-1];
02971 if (it2 <= 0) goto L1250;
02972
02973 L1210:
02974 mnscan();
02975 return;
02976 L1250:
02977 Printf(" PARAMETER %4d NOT VARIABLE.",iext);
02978 ierflg = 3;
02979 return;
02980
02981 L1300:
02982 ke1 = Int_t(fWord7[0]);
02983 ke2 = Int_t(fWord7[1]);
02984 if (ke1 == 0) {
02985 if (fNpar == 2) {
02986 ke1 = fNexofi[0];
02987 ke2 = fNexofi[1];
02988 } else {
02989 Printf("%s: NO PARAMETERS REQUESTED ",(const char*)fCword);
02990 ierflg = 3;
02991 return;
02992 }
02993 }
02994 fNfcnmx = 1000;
02995 mncntr(ke1-1, ke2-1, ierrf);
02996 if (ierrf > 0) ierflg = 3;
02997 return;
02998
02999 L1400:
03000 mnhess();
03001 mnwerr();
03002 if (fISW[4] >= 0) mnprin(2, fAmin);
03003 if (fISW[4] >= 1) mnmatu(1);
03004 return;
03005
03006 L1500:
03007 mnsave();
03008 return;
03009
03010 L1600:
03011 mncuve();
03012 mnimpr();
03013 if (fLnewmn) goto L400;
03014 ierflg = 4;
03015 return;
03016
03017 L1700:
03018 iflag = Int_t(fWord7[0]);
03019 nparx = fNpar;
03020 f = fUndefi;
03021 Eval(nparx, fGin, f, fU, iflag); ++fNfcn;
03022 nowprt = 0;
03023 if (f != fUndefi) {
03024 if (fAmin == fUndefi) {
03025 fAmin = f;
03026 nowprt = 1;
03027 } else if (f < fAmin) {
03028 fAmin = f;
03029 nowprt = 1;
03030 }
03031 if (fISW[4] >= 0 && iflag <= 5 && nowprt == 1) {
03032 mnprin(5, fAmin);
03033 }
03034 if (iflag == 3) fFval3 = f;
03035 }
03036 if (iflag > 5) mnrset(1);
03037 return;
03038
03039 L1800:
03040
03041 return;
03042
03043 L1900:
03044 it = Int_t(fWord7[0]);
03045 if (fFval3 != fAmin && it == 0) {
03046 iflag = 3;
03047 if (fISW[4] >= 0) Printf(" CALL TO USER FUNCTION WITH IFLAG = 3");
03048 nparx = fNpar;
03049 Eval(nparx, fGin, f, fU, iflag); ++fNfcn;
03050 }
03051 ierflg = 11;
03052 if (fCword(0,3) == "END") ierflg = 10;
03053 if (fCword(0,3) == "RET") ierflg = 12;
03054 return;
03055
03056 L2200:
03057 mncler();
03058 if (fISW[4] >= 1) {
03059 Printf(" MINUIT MEMORY CLEARED. NO PARAMETERS NOW DEFINED.");
03060 }
03061 return;
03062
03063 L2300:
03064 kcol = 0;
03065 for (icol = 5; icol <= lk; ++icol) {
03066 if (fCword[icol-1] == ' ') continue;
03067 kcol = icol;
03068 goto L2320;
03069 }
03070 L2320:
03071 if (kcol == 0) comd = "* ";
03072 else comd = fCword(kcol-1,lk-kcol+1);
03073 mnhelp(comd);
03074 return;
03075
03076 L2400:
03077 fEpsi = fUp*.05;
03078 ke1 = Int_t(fWord7[0]);
03079 ke2 = Int_t(fWord7[1]);
03080 if (ke1 == 0 && fNpar == 2) {
03081 ke1 = fNexofi[0];
03082 ke2 = fNexofi[1];
03083 }
03084 nptu = Int_t(fWord7[2]);
03085 if (nptu <= 0) nptu = 20;
03086 if (nptu > 101) nptu = 101;
03087 fNfcnmx = (nptu + 5)*100*(fNpar + 1);
03088 mncont(ke1-1, ke2-1, nptu, xptu, yptu, ierrf);
03089 if (ierrf < nptu) ierflg = 4;
03090 if (ierrf == -1) ierflg = 3;
03091 return;
03092
03093 L2600:
03094 step = fWord7[0];
03095 if (step <= 0) step = 2;
03096 rno = 0;
03097 izero = 0;
03098 for (i = 1; i <= fNpar; ++i) {
03099 mnrn15(rno, izero);
03100 rno = rno*2 - 1;
03101 fX[i-1] += rno*step*fWerr[i-1];
03102 }
03103 mninex(fX);
03104 mnamin();
03105 mnrset(0);
03106 return;
03107
03108 L3300:
03109 Printf(" BLANK COMMAND IGNORED.");
03110 ierflg = 1;
03111 return;
03112
03113
03114 L3400:
03115 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*");
03116 ierflg = 3;
03117 return;
03118
03119 L3500:
03120 cneway = "SET PRInt ";
03121 goto L3100;
03122
03123 L3600:
03124 cneway = "SET GRAd ";
03125 goto L3100;
03126
03127 L3700:
03128 cneway = "SHOW COVar";
03129 goto L3100;
03130
03131 L3800:
03132 cneway = "SET ERRdef";
03133 goto L3100;
03134
03135 L3900:
03136 cneway = "SET LIMits";
03137 goto L3100;
03138
03139 L4000:
03140 cneway = "SAVE ";
03141
03142 L3100:
03143 Printf(" OBSOLETE COMMAND:%s PLEASE USE:%s",(const char*)fCword
03144 ,(const char*)cneway);
03145 fCword = cneway;
03146 if (fCword == "SAVE ") goto L1500;
03147 goto L700;
03148
03149 }
03150
03151
03152 void TMinuit::mnexin(Double_t *pint)
03153 {
03154
03155
03156
03157
03158
03159
03160 Double_t pinti;
03161 Int_t iint, iext;
03162
03163 fLimset = kFALSE;
03164 for (iint = 1; iint <= fNpar; ++iint) {
03165 iext = fNexofi[iint-1];
03166 mnpint(fU[iext-1], iext-1, pinti);
03167 pint[iint-1] = pinti;
03168 }
03169 }
03170
03171
03172 void TMinuit::mnfixp(Int_t iint1, Int_t &ierr)
03173 {
03174
03175
03176
03177
03178
03179
03180 Double_t yyover;
03181 Int_t kold, nold, ndex, knew, iext, i, j, m, n, lc, ik;
03182
03183
03184 ierr = 0;
03185 Int_t iint = iint1+1;
03186 if (iint > fNpar || iint <= 0) {
03187 ierr = 1;
03188 Printf(" MINUIT ERROR. ARGUMENT TO MNFIXP=%4d",iint);
03189 return;
03190 }
03191 iext = fNexofi[iint-1];
03192 if (fNpfix >= fMaxpar) {
03193 ierr = 1;
03194 Printf(" MINUIT CANNOT FIX PARAMETER %4d MAXIMUM NUMBER THAT CAN BE FIXED IS %d",iext,fMaxpar);
03195 return;
03196 }
03197
03198
03199 fNiofex[iext-1] = 0;
03200 nold = fNpar;
03201 --fNpar;
03202
03203
03204 ++fNpfix;
03205 fIpfix[fNpfix-1] = iext;
03206 lc = iint;
03207 fXs[fNpfix-1] = fX[lc-1];
03208 fXts[fNpfix-1] = fXt[lc-1];
03209 fDirins[fNpfix-1] = fWerr[lc-1];
03210 fGrds[fNpfix-1] = fGrd[lc-1];
03211 fG2s[fNpfix-1] = fG2[lc-1];
03212 fGsteps[fNpfix-1] = fGstep[lc-1];
03213
03214 for (ik = iext + 1; ik <= fNu; ++ik) {
03215 if (fNiofex[ik-1] > 0) {
03216 lc = fNiofex[ik-1] - 1;
03217 fNiofex[ik-1] = lc;
03218 fNexofi[lc-1] = ik;
03219 fX[lc-1] = fX[lc];
03220 fXt[lc-1] = fXt[lc];
03221 fDirin[lc-1] = fDirin[lc];
03222 fWerr[lc-1] = fWerr[lc];
03223 fGrd[lc-1] = fGrd[lc];
03224 fG2[lc-1] = fG2[lc];
03225 fGstep[lc-1] = fGstep[lc];
03226 }
03227 }
03228 if (fISW[1] <= 0) return;
03229
03230 if (fNpar <= 0) return;
03231 for (i = 1; i <= nold; ++i) {
03232 m = TMath::Max(i,iint);
03233 n = TMath::Min(i,iint);
03234 ndex = m*(m-1) / 2 + n;
03235 fFIXPyy[i-1] = fVhmat[ndex-1];
03236 }
03237 yyover = 1 / fFIXPyy[iint-1];
03238 knew = 0;
03239 kold = 0;
03240 for (i = 1; i <= nold; ++i) {
03241 for (j = 1; j <= i; ++j) {
03242 ++kold;
03243 if (j == iint || i == iint) continue;
03244 ++knew;
03245 fVhmat[knew-1] = fVhmat[kold-1] - fFIXPyy[j-1]*fFIXPyy[i-1]*yyover;
03246 }
03247 }
03248 }
03249
03250
03251 void TMinuit::mnfree(Int_t k)
03252 {
03253
03254
03255
03256
03257
03258
03259
03260
03261
03262
03263
03264
03265
03266
03267
03268 Double_t grdv, xv, dirinv, g2v, gstepv, xtv;
03269 Int_t i, ipsav, ka, lc, ik, iq, ir, is;
03270
03271 if (k > 1) {
03272 Printf(" CALL TO MNFREE IGNORED. ARGUMENT GREATER THAN ONE");
03273 }
03274 if (fNpfix < 1) {
03275 Printf(" CALL TO MNFREE IGNORED. THERE ARE NO FIXED PARAMETERS");
03276 }
03277 if (k == 1 || k == 0) goto L40;
03278
03279
03280 ka = abs(k);
03281 if (fNiofex[ka-1] == 0) goto L15;
03282 Printf(" IGNORED. PARAMETER SPECIFIED IS ALREADY VARIABLE.");
03283 return;
03284 L15:
03285 if (fNpfix < 1) goto L21;
03286 for (ik = 1; ik <= fNpfix; ++ik) { if (fIpfix[ik-1] == ka) goto L24; }
03287 L21:
03288 Printf(" PARAMETER %4d NOT FIXED. CANNOT BE RELEASED.",ka);
03289 return;
03290 L24:
03291 if (ik == fNpfix) goto L40;
03292
03293
03294 ipsav = ka;
03295 xv = fXs[ik-1];
03296 xtv = fXts[ik-1];
03297 dirinv = fDirins[ik-1];
03298 grdv = fGrds[ik-1];
03299 g2v = fG2s[ik-1];
03300 gstepv = fGsteps[ik-1];
03301 for (i = ik + 1; i <= fNpfix; ++i) {
03302 fIpfix[i-2] = fIpfix[i-1];
03303 fXs[i-2] = fXs[i-1];
03304 fXts[i-2] = fXts[i-1];
03305 fDirins[i-2] = fDirins[i-1];
03306 fGrds[i-2] = fGrds[i-1];
03307 fG2s[i-2] = fG2s[i-1];
03308 fGsteps[i-2] = fGsteps[i-1];
03309 }
03310 fIpfix[fNpfix-1] = ipsav;
03311 fXs[fNpfix-1] = xv;
03312 fXts[fNpfix-1] = xtv;
03313 fDirins[fNpfix-1] = dirinv;
03314 fGrds[fNpfix-1] = grdv;
03315 fG2s[fNpfix-1] = g2v;
03316 fGsteps[fNpfix-1] = gstepv;
03317
03318 L40:
03319 if (fNpfix < 1) goto L300;
03320 ir = fIpfix[fNpfix-1];
03321 is = 0;
03322 for (ik = fNu; ik >= ir; --ik) {
03323 if (fNiofex[ik-1] > 0) {
03324 lc = fNiofex[ik-1] + 1;
03325 is = lc - 1;
03326 fNiofex[ik-1] = lc;
03327 fNexofi[lc-1] = ik;
03328 fX[lc-1] = fX[lc-2];
03329 fXt[lc-1] = fXt[lc-2];
03330 fDirin[lc-1] = fDirin[lc-2];
03331 fWerr[lc-1] = fWerr[lc-2];
03332 fGrd[lc-1] = fGrd[lc-2];
03333 fG2[lc-1] = fG2[lc-2];
03334 fGstep[lc-1] = fGstep[lc-2];
03335 }
03336 }
03337 ++fNpar;
03338 if (is == 0) is = fNpar;
03339 fNiofex[ir-1] = is;
03340 fNexofi[is-1] = ir;
03341 iq = fNpfix;
03342 fX[is-1] = fXs[iq-1];
03343 fXt[is-1] = fXts[iq-1];
03344 fDirin[is-1] = fDirins[iq-1];
03345 fWerr[is-1] = fDirins[iq-1];
03346 fGrd[is-1] = fGrds[iq-1];
03347 fG2[is-1] = fG2s[iq-1];
03348 fGstep[is-1] = fGsteps[iq-1];
03349 --fNpfix;
03350 fISW[1] = 0;
03351 fDcovar = 1;
03352 if (fISW[4] - fItaur >= 1) {
03353 Printf(" PARAMETER %4d %s RESTORED TO VARIABLE.",ir,
03354 (const char*)fCpnam[ir-1]);
03355 }
03356 if (k == 0) goto L40;
03357 L300:
03358
03359 mnexin(fX);
03360 }
03361
03362
03363 void TMinuit::mngrad()
03364 {
03365
03366
03367
03368
03369
03370
03371
03372
03373
03374
03375 Double_t fzero, err;
03376 Int_t i, nparx, lc, istsav;
03377 Bool_t lnone;
03378 static TString cwd = " ";
03379
03380 fISW[2] = 1;
03381 nparx = fNpar;
03382 if (fWord7[0] > 0) goto L2000;
03383
03384
03385 for (i = 1; i <= fNu; ++i) { fGin[i-1] = fUndefi; }
03386 mninex(fX);
03387 Eval(nparx, fGin, fzero, fU, 2); ++fNfcn;
03388 mnderi();
03389 for (i = 1; i <= fNpar; ++i) { fGRADgf[i-1] = fGrd[i-1]; }
03390
03391 fISW[2] = 0;
03392 istsav = fIstrat;
03393 fIstrat = 2;
03394 mnhes1();
03395 fIstrat = istsav;
03396 Printf(" CHECK OF GRADIENT CALCULATION IN FCN");
03397 Printf(" PARAMETER G(IN FCN) G(MINUIT) DG(MINUIT) AGREEMENT");
03398 fISW[2] = 1;
03399 lnone = kFALSE;
03400 for (lc = 1; lc <= fNpar; ++lc) {
03401 i = fNexofi[lc-1];
03402 cwd = "GOOD";
03403 err = fDgrd[lc-1];
03404 if (TMath::Abs(fGRADgf[lc-1] - fGrd[lc-1]) > err) cwd = " BAD";
03405 if (fGin[i-1] == fUndefi) {
03406 cwd = "NONE";
03407 lnone = kTRUE;
03408 fGRADgf[lc-1] = 0;
03409 }
03410 if (cwd != "GOOD") fISW[2] = 0;
03411 Printf(" %5d %10s%12.4e%12.4e%12.4e %s",i
03412 ,(const char*)fCpnam[i-1]
03413 ,fGRADgf[lc-1],fGrd[lc-1],err,(const char*)cwd);
03414 }
03415 if (lnone) {
03416 Printf(" AGREEMENT=NONE MEANS FCN DID NOT CALCULATE THE DERIVATIVE");
03417 }
03418 if (fISW[2] == 0) {
03419 Printf(" MINUIT DOES NOT ACCEPT DERIVATIVE CALCULATIONS BY FCN");
03420 Printf(" TO FORCE ACCEPTANCE, ENTER *SET GRAD 1*");
03421 }
03422
03423 L2000:
03424 return;
03425 }
03426
03427
03428 void TMinuit::mnhelp(const char *command)
03429 {
03430
03431 TString comd = command;
03432 mnhelp(comd);
03433 }
03434
03435
03436 void TMinuit::mnhelp(TString comd)
03437 {
03438
03439
03440
03441
03442
03443
03444
03445
03446
03447
03448
03449
03450
03451
03452
03453
03454
03455
03456 comd.ToUpper();
03457 if( comd.Length() == 0 || comd[0] == '*' || comd[0] == '?' || comd[0] == 0 || comd=="HELP" ) {
03458 Printf(" ==>List of MINUIT Interactive commands:");
03459 Printf(" CLEar Reset all parameter names and values undefined");
03460 Printf(" CONtour Make contour map of the user function");
03461 Printf(" EXIT Exit from Interactive Minuit");
03462 Printf(" FIX Cause parameter(s) to remain constant");
03463 Printf(" HESse Calculate the Hessian or error matrix.");
03464 Printf(" IMPROVE Search for a new minimum around current minimum");
03465 Printf(" MIGrad Minimize by the method of Migrad");
03466 Printf(" MINImize MIGRAD + SIMPLEX method if Migrad fails");
03467 Printf(" MINOs Exact (non-linear) parameter error analysis");
03468 Printf(" MNContour Calculate one MINOS function contour");
03469 Printf(" PARameter Define or redefine new parameters and values");
03470 Printf(" RELease Make previously FIXed parameters variable again");
03471 Printf(" REStore Release last parameter fixed");
03472 Printf(" SAVe Save current parameter values on a file");
03473 Printf(" SCAn Scan the user function by varying parameters");
03474 Printf(" SEEk Minimize by the method of Monte Carlo");
03475 Printf(" SET Set various MINUIT constants or conditions");
03476 Printf(" SHOw Show values of current constants or conditions");
03477 Printf(" SIMplex Minimize by the method of Simplex");
03478 goto L99;
03479 }
03480
03481
03482
03483
03484
03485
03486 if( !strncmp(comd.Data(),"CLE",3) ) {
03487 Printf(" ***>CLEAR");
03488 Printf(" Resets all parameter names and values to undefined.");
03489 Printf(" Must normally be followed by a PARameters command or ");
03490 Printf(" equivalent, in order to define parameter values.");
03491 goto L99;
03492 }
03493
03494
03495
03496
03497
03498 if( !strncmp(comd.Data(),"CON",3) ) {
03499 Printf(" ***>CONTOUR <par1> <par2> [devs] [ngrid]");
03500 Printf(" Instructs Minuit to trace contour lines of the user function");
03501 Printf(" with respect to the two parameters whose external numbers");
03502 Printf(" are <par1> and <par2>.");
03503 Printf(" Other variable parameters of the function, if any, will have");
03504 Printf(" their values fixed at the current values during the contour");
03505 Printf(" tracing. The optional parameter [devs] (default value 2.)");
03506 Printf(" gives the number of standard deviations in each parameter");
03507 Printf(" which should lie entirely within the plotting area.");
03508 Printf(" Optional parameter [ngrid] (default value 25 unless page");
03509 Printf(" size is too small) determines the resolution of the plot,");
03510 Printf(" i.e. the number of rows and columns of the grid at which the");
03511 Printf(" function will be evaluated. [See also MNContour.]");
03512 goto L99;
03513 }
03514
03515
03516
03517
03518
03519 if( !strncmp(comd.Data(),"END",3) ) {
03520 Printf(" ***>END");
03521 Printf(" Signals the end of a data block (i.e., the end of a fit),");
03522 Printf(" and implies that execution should continue, because another");
03523 Printf(" Data Block follows. A Data Block is a set of Minuit data");
03524 Printf(" consisting of");
03525 Printf(" (1) A Title,");
03526 Printf(" (2) One or more Parameter Definitions,");
03527 Printf(" (3) A blank line, and");
03528 Printf(" (4) A set of Minuit Commands.");
03529 Printf(" The END command is used when more than one Data Block is to");
03530 Printf(" be used with the same FCN function. It first causes Minuit");
03531 Printf(" to issue a CALL FCN with IFLAG=3, in order to allow FCN to");
03532 Printf(" perform any calculations associated with the final fitted");
03533 Printf(" parameter values, unless a CALL FCN 3 command has already");
03534 Printf(" been executed at the current FCN value.");
03535 goto L99;
03536 }
03537
03538
03539
03540
03541
03542 if( !strncmp(comd.Data(),"EXI",3) ) {
03543 Printf(" ***>EXIT");
03544 Printf(" Signals the end of execution.");
03545 Printf(" The EXIT command first causes Minuit to issue a CALL FCN");
03546 Printf(" with IFLAG=3, to allow FCN to perform any calculations");
03547 Printf(" associated with the final fitted parameter values, unless a");
03548 Printf(" CALL FCN 3 command has already been executed.");
03549 goto L99;
03550 }
03551
03552
03553
03554
03555
03556 if( !strncmp(comd.Data(),"FIX",3) ) {
03557 Printf(" ***>FIX} <parno> [parno] ... [parno]");
03558 Printf(" Causes parameter(s) <parno> to be removed from the list of");
03559 Printf(" variable parameters, and their value(s) will remain constant");
03560 Printf(" during subsequent minimizations, etc., until another command");
03561 Printf(" changes their value(s) or status.");
03562 goto L99;
03563 }
03564
03565
03566
03567
03568
03569 if( !strncmp(comd.Data(),"HES",3) ) {
03570 Printf(" ***>HESse [maxcalls]");
03571 Printf(" Calculate, by finite differences, the Hessian or error matrix.");
03572 Printf(" That is, it calculates the full matrix of second derivatives");
03573 Printf(" of the function with respect to the currently variable");
03574 Printf(" parameters, and inverts it, printing out the resulting error");
03575 Printf(" matrix. The optional argument [maxcalls] specifies the");
03576 Printf(" (approximate) maximum number of function calls after which");
03577 Printf(" the calculation will be stopped.");
03578 goto L99;
03579 }
03580
03581
03582
03583
03584
03585 if( !strncmp(comd.Data(),"IMP",3) ) {
03586 Printf(" ***>IMPROVE [maxcalls]");
03587 Printf(" If a previous minimization has converged, and the current");
03588 Printf(" values of the parameters therefore correspond to a local");
03589 Printf(" minimum of the function, this command requests a search for");
03590 Printf(" additional distinct local minima.");
03591 Printf(" The optional argument [maxcalls] specifies the (approximate");
03592 Printf(" maximum number of function calls after which the calculation");
03593 Printf(" will be stopped.");
03594 goto L99;
03595 }
03596
03597
03598
03599
03600
03601 if( !strncmp(comd.Data(),"MIG",3) ) {
03602 Printf(" ***>MIGrad [maxcalls] [tolerance]");
03603 Printf(" Causes minimization of the function by the method of Migrad,");
03604 Printf(" the most efficient and complete single method, recommended");
03605 Printf(" for general functions (see also MINImize).");
03606 Printf(" The minimization produces as a by-product the error matrix");
03607 Printf(" of the parameters, which is usually reliable unless warning");
03608 Printf(" messages are produced.");
03609 Printf(" The optional argument [maxcalls] specifies the (approximate)");
03610 Printf(" maximum number of function calls after which the calculation");
03611 Printf(" will be stopped even if it has not yet converged.");
03612 Printf(" The optional argument [tolerance] specifies required tolerance");
03613 Printf(" on the function value at the minimum.");
03614 Printf(" The default tolerance is 0.1, and the minimization will stop");
03615 Printf(" when the estimated vertical distance to the minimum (EDM) is");
03616 Printf(" less than 0.001*[tolerance]*UP (see [SET ERRordef]).");
03617 goto L99;
03618 }
03619
03620
03621
03622
03623
03624 if( !strncmp(comd.Data(),"MINI",4) ) {
03625 Printf(" ***>MINImize [maxcalls] [tolerance]");
03626 Printf(" Causes minimization of the function by the method of Migrad,");
03627 Printf(" as does the MIGrad command, but switches to the SIMplex method");
03628 Printf(" if Migrad fails to converge. Arguments are as for MIGrad.");
03629 Printf(" Note that command requires four characters to be unambiguous.");
03630 goto L99;
03631 }
03632
03633
03634
03635
03636
03637 if( !strncmp(comd.Data(),"MIN0",4) ) {
03638 Printf(" ***>MINOs [maxcalls] [parno] [parno] ...");
03639 Printf(" Causes a Minos error analysis to be performed on the parameters");
03640 Printf(" whose numbers [parno] are specified. If none are specified,");
03641 Printf(" Minos errors are calculated for all variable parameters.");
03642 Printf(" Minos errors may be expensive to calculate, but are very");
03643 Printf(" reliable since they take account of non-linearities in the");
03644 Printf(" problem as well as parameter correlations, and are in general");
03645 Printf(" asymmetric.");
03646 Printf(" The optional argument [maxcalls] specifies the (approximate)");
03647 Printf(" maximum number of function calls per parameter requested,");
03648 Printf(" after which the calculation will stop for that parameter.");
03649 goto L99;
03650 }
03651
03652
03653
03654
03655
03656 if( !strncmp(comd.Data(),"MNC",3) ) {
03657 Printf(" ***>MNContour <par1> <par2> [npts]");
03658 Printf(" Calculates one function contour of FCN with respect to");
03659 Printf(" parameters par1 and par2, with FCN minimized always with");
03660 Printf(" respect to all other NPAR-2 variable parameters (if any).");
03661 Printf(" Minuit will try to find npts points on the contour (default 20)");
03662 Printf(" If only two parameters are variable at the time, it is not");
03663 Printf(" necessary to specify their numbers. To calculate more than");
03664 Printf(" one contour, it is necessary to SET ERRordef to the appropriate");
03665 Printf(" value and issue the MNContour command for each contour.");
03666 goto L99;
03667 }
03668
03669
03670
03671
03672
03673 if( !strncmp(comd.Data(),"PAR",3) ) {
03674 Printf(" ***>PARameters");
03675 Printf(" followed by one or more parameter definitions.");
03676 Printf(" Parameter definitions are of the form:");
03677 Printf(" <number> ''name'' <value> <step> [lolim] [uplim] ");
03678 Printf(" for example:");
03679 Printf(" 3 ''K width'' 1.2 0.1");
03680 Printf(" the last definition is followed by a blank line or a zero.");
03681 goto L99;
03682 }
03683
03684
03685
03686
03687
03688 if( !strncmp(comd.Data(),"REL",3) ) {
03689 Printf(" ***>RELease <parno> [parno] ... [parno]");
03690 Printf(" If <parno> is the number of a previously variable parameter");
03691 Printf(" which has been fixed by a command: FIX <parno>, then that");
03692 Printf(" parameter will return to variable status. Otherwise a warning");
03693 Printf(" message is printed and the command is ignored.");
03694 Printf(" Note that this command operates only on parameters which were");
03695 Printf(" at one time variable and have been FIXed. It cannot make");
03696 Printf(" constant parameters variable; that must be done by redefining");
03697 Printf(" the parameter with a PARameters command.");
03698 goto L99;
03699 }
03700
03701
03702
03703
03704
03705 if( !strncmp(comd.Data(),"RES",3) ) {
03706 Printf(" ***>REStore [code]");
03707 Printf(" If no [code] is specified, this command restores all previously");
03708 Printf(" FIXed parameters to variable status. If [code]=1, then only");
03709 Printf(" the last parameter FIXed is restored to variable status.");
03710 Printf(" If code is neither zero nor one, the command is ignored.");
03711 goto L99;
03712 }
03713
03714
03715
03716
03717
03718 if( !strncmp(comd.Data(),"RET",3) ) {
03719 Printf(" ***>RETURN");
03720 Printf(" Signals the end of a data block, and instructs Minuit to return");
03721 Printf(" to the program which called it. The RETurn command first");
03722 Printf(" causes Minuit to CALL FCN with IFLAG=3, in order to allow FCN");
03723 Printf(" to perform any calculations associated with the final fitted");
03724 Printf(" parameter values, unless a CALL FCN 3 command has already been");
03725 Printf(" executed at the current FCN value.");
03726 goto L99;
03727 }
03728
03729
03730
03731
03732
03733 if( !strncmp(comd.Data(),"SAV",3) ) {
03734 Printf(" ***>SAVe");
03735 Printf(" Causes the current parameter values to be saved on a file in");
03736 Printf(" such a format that they can be read in again as Minuit");
03737 Printf(" parameter definitions. If the covariance matrix exists, it is");
03738 Printf(" also output in such a format. The unit number is by default 7,");
03739 Printf(" or that specified by the user in his call to MINTIO or");
03740 Printf(" MNINIT. The user is responsible for opening the file previous");
03741 Printf(" to issuing the [SAVe] command (except where this can be done");
03742 Printf(" interactively).");
03743 goto L99;
03744 }
03745
03746
03747
03748
03749
03750 if( !strncmp(comd.Data(),"SCA",3) ) {
03751 Printf(" ***>SCAn [parno] [numpts] [from] [to]");
03752 Printf(" Scans the value of the user function by varying parameter");
03753 Printf(" number [parno], leaving all other parameters fixed at the");
03754 Printf(" current value. If [parno] is not specified, all variable");
03755 Printf(" parameters are scanned in sequence.");
03756 Printf(" The number of points [numpts] in the scan is 40 by default,");
03757 Printf(" and cannot exceed 100. The range of the scan is by default");
03758 Printf(" 2 standard deviations on each side of the current best value,");
03759 Printf(" but can be specified as from [from] to [to].");
03760 Printf(" After each scan, if a new minimum is found, the best parameter");
03761 Printf(" values are retained as start values for future scans or");
03762 Printf(" minimizations. The curve resulting from each scan is plotted");
03763 Printf(" on the output unit in order to show the approximate behaviour");
03764 Printf(" of the function.");
03765 Printf(" This command is not intended for minimization, but is sometimes");
03766 Printf(" useful for debugging the user function or finding a");
03767 Printf(" reasonable starting point.");
03768 goto L99;
03769 }
03770
03771
03772
03773
03774
03775 if( !strncmp(comd.Data(),"SEE",3) ) {
03776 Printf(" ***>SEEk [maxcalls] [devs]");
03777 Printf(" Causes a Monte Carlo minimization of the function, by choosing");
03778 Printf(" random values of the variable parameters, chosen uniformly");
03779 Printf(" over a hypercube centered at the current best value.");
03780 Printf(" The region size is by default 3 standard deviations on each");
03781 Printf(" side, but can be changed by specifying the value of [devs].");
03782 goto L99;
03783 }
03784
03785
03786
03787
03788
03789 if( !strncmp(comd.Data(),"SET",3) ) {
03790 Printf(" ***>SET <option_name>");
03791 Printf(" SET BATch");
03792 Printf(" Informs Minuit that it is running in batch mode.");
03793
03794 Printf(" ");
03795 Printf(" SET EPSmachine <accuracy>");
03796 Printf(" Informs Minuit that the relative floating point arithmetic");
03797 Printf(" precision is <accuracy>. Minuit determines the nominal");
03798 Printf(" precision itself, but the SET EPSmachine command can be");
03799 Printf(" used to override Minuit own determination, when the user");
03800 Printf(" knows that the FCN function value is not calculated to");
03801 Printf(" the nominal machine accuracy. Typical values of <accuracy>");
03802 Printf(" are between 10**-5 and 10**-14.");
03803
03804 Printf(" ");
03805 Printf(" SET ERRordef <up>");
03806 Printf(" Sets the value of UP (default value= 1.), defining parameter");
03807 Printf(" errors. Minuit defines parameter errors as the change");
03808 Printf(" in parameter value required to change the function value");
03809 Printf(" by UP. Normally, for chisquared fits UP=1, and for negative");
03810 Printf(" log likelihood, UP=0.5.");
03811
03812 Printf(" ");
03813 Printf(" SET GRAdient [force]");
03814 Printf(" Informs Minuit that the user function is prepared to");
03815 Printf(" calculate its own first derivatives and return their values");
03816 Printf(" in the array GRAD when IFLAG=2 (see specs of FCN).");
03817 Printf(" If [force] is not specified, Minuit will calculate");
03818 Printf(" the FCN derivatives by finite differences at the current");
03819 Printf(" point and compare with the user calculation at that point,");
03820 Printf(" accepting the user values only if they agree.");
03821 Printf(" If [force]=1, Minuit does not do its own derivative");
03822 Printf(" calculation, and uses the derivatives calculated in FCN.");
03823
03824 Printf(" ");
03825 Printf(" SET INPut [unitno] [filename]");
03826 Printf(" Causes Minuit, in data-driven mode only, to read subsequent");
03827 Printf(" commands (or parameter definitions) from a different input");
03828 Printf(" file. If no [unitno] is specified, reading reverts to the");
03829 Printf(" previous input file, assuming that there was one.");
03830 Printf(" If [unitno] is specified, and that unit has not been opened,");
03831 Printf(" then Minuit attempts to open the file [filename]} if a");
03832 Printf(" name is specified. If running in interactive mode and");
03833 Printf(" [filename] is not specified and [unitno] is not opened,");
03834 Printf(" Minuit prompts the user to enter a file name.");
03835 Printf(" If the word REWIND is added to the command (note:no blanks");
03836 Printf(" between INPUT and REWIND), the file is rewound before");
03837 Printf(" reading. Note that this command is implemented in standard");
03838 Printf(" Fortran 77 and the results may depend on the system;");
03839 Printf(" for example, if a filename is given under VM/CMS, it must");
03840 Printf(" be preceeded by a slash.");
03841
03842 Printf(" ");
03843 Printf(" SET INTeractive");
03844 Printf(" Informs Minuit that it is running interactively.");
03845
03846 Printf(" ");
03847 Printf(" SET LIMits [parno] [lolim] [uplim]");
03848 Printf(" Allows the user to change the limits on one or all");
03849 Printf(" parameters. If no arguments are specified, all limits are");
03850 Printf(" removed from all parameters. If [parno] alone is specified,");
03851 Printf(" limits are removed from parameter [parno].");
03852 Printf(" If all arguments are specified, then parameter [parno] will");
03853 Printf(" be bounded between [lolim] and [uplim].");
03854 Printf(" Limits can be specified in either order, Minuit will take");
03855 Printf(" the smaller as [lolim] and the larger as [uplim].");
03856 Printf(" However, if [lolim] is equal to [uplim], an error condition");
03857 Printf(" results.");
03858
03859 Printf(" ");
03860 Printf(" SET LINesperpage");
03861 Printf(" Sets the number of lines for one page of output.");
03862 Printf(" Default value is 24 for interactive mode");
03863
03864 Printf(" ");
03865 Printf(" SET NOGradient");
03866 Printf(" The inverse of SET GRAdient, instructs Minuit not to");
03867 Printf(" use the first derivatives calculated by the user in FCN.");
03868
03869 Printf(" ");
03870 Printf(" SET NOWarnings");
03871 Printf(" Supresses Minuit warning messages.");
03872
03873 Printf(" ");
03874 Printf(" SET OUTputfile <unitno>");
03875 Printf(" Instructs Minuit to write further output to unit <unitno>.");
03876
03877 Printf(" ");
03878 Printf(" SET PAGethrow <integer>");
03879 Printf(" Sets the carriage control character for ``new page'' to");
03880 Printf(" <integer>. Thus the value 1 produces a new page, and 0");
03881 Printf(" produces a blank line, on some devices (see TOPofpage)");
03882
03883
03884 Printf(" ");
03885 Printf(" SET PARameter <parno> <value>");
03886 Printf(" Sets the value of parameter <parno> to <value>.");
03887 Printf(" The parameter in question may be variable, fixed, or");
03888 Printf(" constant, but must be defined.");
03889
03890 Printf(" ");
03891 Printf(" SET PRIntout <level>");
03892 Printf(" Sets the print level, determining how much output will be");
03893 Printf(" produced. Allowed values and their meanings are displayed");
03894 Printf(" after a SHOw PRInt command, and are currently <level>=:");
03895 Printf(" [-1] no output except from SHOW commands");
03896 Printf(" [0] minimum output");
03897 Printf(" [1] default value, normal output");
03898 Printf(" [2] additional output giving intermediate results.");
03899 Printf(" [3] maximum output, showing progress of minimizations.");
03900 Printf(" Note: See also the SET WARnings command.");
03901
03902 Printf(" ");
03903 Printf(" SET RANdomgenerator <seed>");
03904 Printf(" Sets the seed of the random number generator used in SEEk.");
03905 Printf(" This can be any integer between 10000 and 900000000, for");
03906 Printf(" example one which was output from a SHOw RANdom command of");
03907 Printf(" a previous run.");
03908
03909 Printf(" ");
03910 Printf(" SET STRategy <level>");
03911 Printf(" Sets the strategy to be used in calculating first and second");
03912 Printf(" derivatives and in certain minimization methods.");
03913 Printf(" In general, low values of <level> mean fewer function calls");
03914 Printf(" and high values mean more reliable minimization.");
03915 Printf(" Currently allowed values are 0, 1 (default), and 2.");
03916
03917 Printf(" ");
03918 Printf(" SET TITle");
03919 Printf(" Informs Minuit that the next input line is to be considered");
03920 Printf(" the (new) title for this task or sub-task. This is for");
03921 Printf(" the convenience of the user in reading his output.");
03922
03923 Printf(" ");
03924 Printf(" SET WARnings");
03925 Printf(" Instructs Minuit to output warning messages when suspicious");
03926 Printf(" conditions arise which may indicate unreliable results.");
03927 Printf(" This is the default.");
03928
03929 Printf(" ");
03930 Printf(" SET WIDthpage");
03931 Printf(" Informs Minuit of the output page width.");
03932 Printf(" Default values are 80 for interactive jobs");
03933 goto L99;
03934 }
03935
03936
03937
03938
03939
03940 if( !strncmp(comd.Data(),"SHO",3) ) {
03941 Printf(" ***>SHOw <option_name>");
03942 Printf(" All SET XXXX commands have a corresponding SHOw XXXX command.");
03943 Printf(" In addition, the SHOw commands listed starting here have no");
03944 Printf(" corresponding SET command for obvious reasons.");
03945
03946 Printf(" ");
03947 Printf(" SHOw CORrelations");
03948 Printf(" Calculates and prints the parameter correlations from the");
03949 Printf(" error matrix.");
03950
03951 Printf(" ");
03952 Printf(" SHOw COVariance");
03953 Printf(" Prints the (external) covariance (error) matrix.");
03954
03955 Printf(" ");
03956 Printf(" SHOw EIGenvalues");
03957 Printf(" Calculates and prints the eigenvalues of the covariance");
03958 Printf(" matrix.");
03959
03960 Printf(" ");
03961 Printf(" SHOw FCNvalue");
03962 Printf(" Prints the current value of FCN.");
03963 goto L99;
03964 }
03965
03966
03967
03968
03969
03970 if( !strncmp(comd.Data(),"SIM",3) ) {
03971 Printf(" ***>SIMplex [maxcalls] [tolerance]");
03972 Printf(" Performs a function minimization using the simplex method of");
03973 Printf(" Nelder and Mead. Minimization terminates either when the");
03974 Printf(" function has been called (approximately) [maxcalls] times,");
03975 Printf(" or when the estimated vertical distance to minimum (EDM) is");
03976 Printf(" less than [tolerance].");
03977 Printf(" The default value of [tolerance] is 0.1*UP(see SET ERRordef).");
03978 goto L99;
03979 }
03980
03981
03982
03983
03984
03985 if( !strncmp(comd.Data(),"STA",3) ) {
03986 Printf(" ***>STAndard");
03987 goto L99;
03988 }
03989
03990
03991
03992
03993
03994 if( !strncmp(comd.Data(),"STO",3) ) {
03995 Printf(" ***>STOP");
03996 Printf(" Same as EXIT.");
03997 goto L99;
03998 }
03999
04000
04001
04002
04003
04004 if( !strncmp(comd.Data(),"TOP",3) ) {
04005 Printf(" ***>TOPofpage");
04006 Printf(" Causes Minuit to write the character specified in a");
04007 Printf(" SET PAGethrow command (default = 1) to column 1 of the output");
04008 Printf(" file, which may or may not position your output medium to");
04009 Printf(" the top of a page depending on the device and system.");
04010 goto L99;
04011 }
04012
04013 Printf(" Unknown MINUIT command. Type HELP for list of commands.");
04014
04015 L99:
04016 return;
04017 }
04018
04019
04020 void TMinuit::mnhess()
04021 {
04022
04023
04024
04025
04026
04027
04028
04029
04030
04031
04032 Double_t dmin_, dxdi, elem, wint, tlrg2, d, dlast, ztemp, g2bfor;
04033 Double_t df, aimsag, fs1, tlrstp, fs2, stpinm, g2i, sag=0, xtf, xti, xtj;
04034 Int_t icyc, ncyc, ndex, idrv, iext, npar2, i, j, ifail, npard, nparx, id, multpy;
04035 Bool_t ldebug;
04036
04037 ldebug = fIdbg[3] >= 1;
04038 if (fAmin == fUndefi) {
04039 mnamin();
04040 }
04041 if (fIstrat <= 0) {
04042 ncyc = 3;
04043 tlrstp = .5;
04044 tlrg2 = .1;
04045 } else if (fIstrat == 1) {
04046 ncyc = 5;
04047 tlrstp = .3;
04048 tlrg2 = .05;
04049 } else {
04050 ncyc = 7;
04051 tlrstp = .1;
04052 tlrg2 = .02;
04053 }
04054 if (fISW[4] >= 2 || ldebug) {
04055 Printf(" START COVARIANCE MATRIX CALCULATION.");
04056 }
04057 fCfrom = "HESSE ";
04058 fNfcnfr = fNfcn;
04059 fCstatu = "OK ";
04060 npard = fNpar;
04061
04062 mninex(fX);
04063 nparx = fNpar;
04064 Eval(nparx, fGin, fs1, fU, 4); ++fNfcn;
04065 if (fs1 != fAmin) {
04066 df = fAmin - fs1;
04067 mnwarn("D", "MNHESS", Form("function value differs from AMIN by %g",df));
04068 }
04069 fAmin = fs1;
04070 if (ldebug) {
04071 Printf(" PAR D GSTEP D G2 GRD SAG ");
04072 }
04073
04074
04075
04076
04077
04078 aimsag = TMath::Sqrt(fEpsma2)*(TMath::Abs(fAmin) + fUp);
04079
04080 npar2 = fNpar*(fNpar + 1) / 2;
04081 for (i = 1; i <= npar2; ++i) { fVhmat[i-1] = 0; }
04082
04083
04084 idrv = 2;
04085 for (id = 1; id <= npard; ++id) {
04086 i = id + fNpar - npard;
04087 iext = fNexofi[i-1];
04088 if (fG2[i-1] == 0) {
04089 mnwarn("W", "HESSE", Form("Second derivative enters zero, param %d",iext));
04090 wint = fWerr[i-1];
04091 if (fNvarl[iext-1] > 1) {
04092 mndxdi(fX[i-1], i-1, dxdi);
04093 if (TMath::Abs(dxdi) < .001) wint = .01;
04094 else wint /= TMath::Abs(dxdi);
04095 }
04096 fG2[i-1] = fUp / (wint*wint);
04097 }
04098 xtf = fX[i-1];
04099 dmin_ = fEpsma2*8*TMath::Abs(xtf);
04100
04101
04102 d = TMath::Abs(fGstep[i-1]);
04103 int skip50 = 0;
04104 for (icyc = 1; icyc <= ncyc; ++icyc) {
04105
04106 for (multpy = 1; multpy <= 5; ++multpy) {
04107
04108 fX[i-1] = xtf + d;
04109 mninex(fX);
04110 nparx = fNpar;
04111 Eval(nparx, fGin, fs1, fU, 4); ++fNfcn;
04112 fX[i-1] = xtf - d;
04113 mninex(fX);
04114 Eval(nparx, fGin, fs2, fU, 4); ++fNfcn;
04115 fX[i-1] = xtf;
04116 sag = (fs1 + fs2 - fAmin*2)*.5;
04117 if (sag != 0) goto L30;
04118 if (fGstep[i-1] < 0) {
04119 if (d >= .5) goto L26;
04120 d *= 10;
04121 if (d > .5) d = .51;
04122 continue;
04123 }
04124 d *= 10;
04125 }
04126 L26:
04127 mnwarn("W", "HESSE", Form("Second derivative zero for parameter%d",iext));
04128 goto L390;
04129
04130 L30:
04131 g2bfor = fG2[i-1];
04132 fG2[i-1] = sag*2 / (d*d);
04133 fGrd[i-1] = (fs1 - fs2) / (d*2);
04134 if (ldebug) {
04135 Printf("%4d%2d%12.5g%12.5g%12.5g%12.5g",i,idrv,fGstep[i-1],fG2[i-1],fGrd[i-1],sag);
04136 }
04137 if (fGstep[i-1] > 0) fGstep[i-1] = TMath::Abs(d);
04138 else fGstep[i-1] = -TMath::Abs(d);
04139 fDirin[i-1] = d;
04140 fHESSyy[i-1]= fs1;
04141 dlast = d;
04142 d = TMath::Sqrt(aimsag*2 / TMath::Abs(fG2[i-1]));
04143
04144 stpinm = .5;
04145 if (fGstep[i-1] < 0) d = TMath::Min(d,stpinm);
04146 if (d < dmin_) d = dmin_;
04147
04148 if (TMath::Abs((d - dlast) / d) < tlrstp ||
04149 TMath::Abs((fG2[i-1] - g2bfor) / fG2[i-1]) < tlrg2) {
04150 skip50 = 1;
04151 break;
04152 }
04153 d = TMath::Min(d,dlast*102);
04154 d = TMath::Max(d,dlast*.1);
04155 }
04156
04157 if (!skip50)
04158 mnwarn("D", "MNHESS", Form("Second Deriv. SAG,AIM= %d%g%g",iext,sag,aimsag));
04159
04160 ndex = i*(i + 1) / 2;
04161 fVhmat[ndex-1] = fG2[i-1];
04162 }
04163
04164 mninex(fX);
04165
04166 if (fIstrat > 0) mnhes1();
04167 fISW[1] = 3;
04168 fDcovar = 0;
04169
04170
04171 if (fNpar == 1) goto L214;
04172 for (i = 1; i <= fNpar; ++i) {
04173 for (j = 1; j <= i-1; ++j) {
04174 xti = fX[i-1];
04175 xtj = fX[j-1];
04176 fX[i-1] = xti + fDirin[i-1];
04177 fX[j-1] = xtj + fDirin[j-1];
04178 mninex(fX);
04179 Eval(nparx, fGin, fs1, fU, 4); ++fNfcn;
04180 fX[i-1] = xti;
04181 fX[j-1] = xtj;
04182 elem = (fs1 + fAmin - fHESSyy[i-1] - fHESSyy[j-1]) / (
04183 fDirin[i-1]*fDirin[j-1]);
04184 ndex = i*(i-1) / 2 + j;
04185 fVhmat[ndex-1] = elem;
04186 }
04187 }
04188 L214:
04189 mninex(fX);
04190
04191 mnpsdf();
04192 for (i = 1; i <= fNpar; ++i) {
04193 for (j = 1; j <= i; ++j) {
04194 ndex = i*(i-1) / 2 + j;
04195 fP[i + j*fMaxpar - fMaxpar-1] = fVhmat[ndex-1];
04196 fP[j + i*fMaxpar - fMaxpar-1] = fP[i + j*fMaxpar - fMaxpar-1];
04197 }
04198 }
04199 mnvert(fP, fMaxint, fMaxint, fNpar, ifail);
04200 if (ifail > 0) {
04201 mnwarn("W", "HESSE", "Matrix inversion fails.");
04202 goto L390;
04203 }
04204
04205 fEDM = 0;
04206
04207 for (i = 1; i <= fNpar; ++i) {
04208
04209 ndex = i*(i-1) / 2;
04210 for (j = 1; j <= i-1; ++j) {
04211 ++ndex;
04212 ztemp = fP[i + j*fMaxpar - fMaxpar-1]*2;
04213 fEDM += fGrd[i-1]*ztemp*fGrd[j-1];
04214 fVhmat[ndex-1] = ztemp;
04215 }
04216
04217 ++ndex;
04218 fVhmat[ndex-1] = fP[i + i*fMaxpar - fMaxpar-1]*2;
04219 fEDM += fP[i + i*fMaxpar - fMaxpar-1]*(fGrd[i-1]*fGrd[i-1]);
04220 }
04221 if (fISW[4] >= 1 && fISW[1] == 3 && fItaur == 0) {
04222 Printf(" COVARIANCE MATRIX CALCULATED SUCCESSFULLY");
04223 }
04224 goto L900;
04225
04226 L390:
04227 fISW[1] = 1;
04228 fDcovar = 1;
04229 fCstatu = "FAILED ";
04230 if (fISW[4] >= 0) {
04231 Printf(" MNHESS FAILS AND WILL RETURN DIAGONAL MATRIX. ");
04232 }
04233 for (i = 1; i <= fNpar; ++i) {
04234 ndex = i*(i-1) / 2;
04235 for (j = 1; j <= i-1; ++j) {
04236 ++ndex;
04237 fVhmat[ndex-1] = 0;
04238 }
04239 ++ndex;
04240 g2i = fG2[i-1];
04241 if (g2i <= 0) g2i = 1;
04242 fVhmat[ndex-1] = 2 / g2i;
04243 }
04244 L900:
04245 return;
04246 }
04247
04248
04249 void TMinuit::mnhes1()
04250 {
04251
04252
04253
04254
04255
04256
04257
04258 Double_t dmin_, d, dfmin, dgmin=0, change, chgold, grdold=0, epspri;
04259 Double_t fs1, optstp, fs2, grdnew=0, sag, xtf;
04260 Int_t icyc, ncyc=0, idrv, i, nparx;
04261 Bool_t ldebug;
04262
04263 ldebug = fIdbg[5] >= 1;
04264 if (fIstrat <= 0) ncyc = 1;
04265 if (fIstrat == 1) ncyc = 2;
04266 if (fIstrat > 1) ncyc = 6;
04267 idrv = 1;
04268 nparx = fNpar;
04269 dfmin = fEpsma2*4*(TMath::Abs(fAmin) + fUp);
04270
04271 for (i = 1; i <= fNpar; ++i) {
04272 xtf = fX[i-1];
04273 dmin_ = fEpsma2*4*TMath::Abs(xtf);
04274 epspri = fEpsma2 + TMath::Abs(fGrd[i-1]*fEpsma2);
04275 optstp = TMath::Sqrt(dfmin / (TMath::Abs(fG2[i-1]) + epspri));
04276 d = TMath::Abs(fGstep[i-1])*.2;
04277 if (d > optstp) d = optstp;
04278 if (d < dmin_) d = dmin_;
04279 chgold = 1e4;
04280
04281 for (icyc = 1; icyc <= ncyc; ++icyc) {
04282 fX[i-1] = xtf + d;
04283 mninex(fX);
04284 Eval(nparx, fGin, fs1, fU, 4); ++fNfcn;
04285 fX[i-1] = xtf - d;
04286 mninex(fX);
04287 Eval(nparx, fGin, fs2, fU, 4); ++fNfcn;
04288 fX[i-1] = xtf;
04289
04290 sag = (fs1 + fs2 - fAmin*2)*.5;
04291 grdold = fGrd[i-1];
04292 grdnew = (fs1 - fs2) / (d*2);
04293 dgmin = fEpsmac*(TMath::Abs(fs1) + TMath::Abs(fs2)) / d;
04294 if (ldebug) {
04295 Printf("%4d%2d%12.5g%12.5g%12.5g%12.5g%12.5g",i,idrv,fGstep[i-1],d,fG2[i-1],grdnew,sag);
04296 }
04297 if (grdnew == 0) goto L60;
04298 change = TMath::Abs((grdold - grdnew) / grdnew);
04299 if (change > chgold && icyc > 1) goto L60;
04300 chgold = change;
04301 fGrd[i-1] = grdnew;
04302 if (fGstep[i-1] > 0) fGstep[i-1] = TMath::Abs(d);
04303 else fGstep[i-1] = -TMath::Abs(d);
04304
04305 if (change < .05) goto L60;
04306 if (TMath::Abs(grdold - grdnew) < dgmin) goto L60;
04307 if (d < dmin_) {
04308 mnwarn("D", "MNHES1", "Step size too small for 1st drv.");
04309 goto L60;
04310 }
04311 d *= .2;
04312 }
04313
04314 mnwarn("D", "MNHES1", Form("Too many iterations on D1.%g%g",grdold,grdnew));
04315 L60:
04316 fDgrd[i-1] = TMath::Max(dgmin,TMath::Abs(grdold - grdnew));
04317 }
04318
04319 mninex(fX);
04320 }
04321
04322
04323 void TMinuit::mnimpr()
04324 {
04325
04326
04327
04328
04329
04330
04331
04332
04333
04334
04335
04336 static Double_t rnum = 0;
04337
04338
04339 Double_t amax, ycalf, ystar, ystst;
04340 Double_t pb, ep, wg, xi, sigsav, reg, sig2;
04341 Int_t npfn, ndex, loop=0, i, j, ifail, iseed=0;
04342 Int_t jhold, nloop, nparx, nparp1, jh, jl, iswtr;
04343
04344 if (fNpar <= 0) return;
04345 if (fAmin == fUndefi) mnamin();
04346 fCstatu = "UNCHANGED ";
04347 fItaur = 1;
04348 fEpsi = fUp*.1;
04349 npfn = fNfcn;
04350 nloop = Int_t(fWord7[1]);
04351 if (nloop <= 0) nloop = fNpar + 4;
04352 nparx = fNpar;
04353 nparp1 = fNpar + 1;
04354 wg = 1 / Double_t(fNpar);
04355 sigsav = fEDM;
04356 fApsi = fAmin;
04357 iswtr = fISW[4] - 2*fItaur;
04358 for (i = 1; i <= fNpar; ++i) {
04359 fXt[i-1] = fX[i-1];
04360 fIMPRdsav[i-1] = fWerr[i-1];
04361 for (j = 1; j <= i; ++j) {
04362 ndex = i*(i-1) / 2 + j;
04363 fP[i + j*fMaxpar - fMaxpar-1] = fVhmat[ndex-1];
04364 fP[j + i*fMaxpar - fMaxpar-1] = fP[i + j*fMaxpar - fMaxpar-1];
04365 }
04366 }
04367 mnvert(fP, fMaxint, fMaxint, fNpar, ifail);
04368 if (ifail >= 1) goto L280;
04369
04370 for (i = 1; i <= fNpar; ++i) {
04371 ndex = i*(i-1) / 2;
04372 for (j = 1; j <= i; ++j) {
04373 ++ndex;
04374 fVthmat[ndex-1] = fP[i + j*fMaxpar - fMaxpar-1];
04375 }
04376 }
04377 loop = 0;
04378
04379 L20:
04380 for (i = 1; i <= fNpar; ++i) {
04381 fDirin[i-1] = fIMPRdsav[i-1]*2;
04382 mnrn15(rnum, iseed);
04383 fX[i-1] = fXt[i-1] + fDirin[i-1]*2*(rnum - .5);
04384 }
04385 ++loop;
04386 reg = 2;
04387 if (fISW[4] >= 0) {
04388 Printf("START ATTEMPT NO.%2d TO FIND NEW MINIMUM",loop);
04389 }
04390 L30:
04391 mncalf(fX, ycalf);
04392 fAmin = ycalf;
04393
04394 jl = nparp1;
04395 jh = nparp1;
04396 fIMPRy[nparp1-1] = fAmin;
04397 amax = fAmin;
04398 for (i = 1; i <= fNpar; ++i) {
04399 xi = fX[i-1];
04400 mnrn15(rnum, iseed);
04401 fX[i-1] = xi - fDirin[i-1]*(rnum - .5);
04402 mncalf(fX, ycalf);
04403 fIMPRy[i-1] = ycalf;
04404 if (fIMPRy[i-1] < fAmin) {
04405 fAmin = fIMPRy[i-1];
04406 jl = i;
04407 } else if (fIMPRy[i-1] > amax) {
04408 amax = fIMPRy[i-1];
04409 jh = i;
04410 }
04411 for (j = 1; j <= fNpar; ++j) { fP[j + i*fMaxpar - fMaxpar-1] = fX[j-1]; }
04412 fP[i + nparp1*fMaxpar - fMaxpar-1] = xi;
04413 fX[i-1] = xi;
04414 }
04415
04416 fEDM = fAmin;
04417 sig2 = fEDM;
04418
04419 L50:
04420 if (fAmin < 0) goto L95;
04421 if (fISW[1] <= 2) goto L280;
04422 ep = fAmin*.1;
04423 if (sig2 < ep && fEDM < ep) goto L100;
04424 sig2 = fEDM;
04425 if (fNfcn - npfn > fNfcnmx) goto L300;
04426
04427 for (i = 1; i <= fNpar; ++i) {
04428 pb = 0;
04429 for (j = 1; j <= nparp1; ++j) { pb += wg*fP[i + j*fMaxpar - fMaxpar-1]; }
04430 fPbar[i-1] = pb - wg*fP[i + jh*fMaxpar - fMaxpar-1];
04431 fPstar[i-1] = fPbar[i-1]*2 - fP[i + jh*fMaxpar - fMaxpar-1]*1;
04432 }
04433 mncalf(fPstar, ycalf);
04434 ystar = ycalf;
04435 if (ystar >= fAmin) goto L70;
04436
04437 for (i = 1; i <= fNpar; ++i) {
04438 fPstst[i-1] = fPstar[i-1]*2 + fPbar[i- 1]*-1;
04439 }
04440 mncalf(fPstst, ycalf);
04441 ystst = ycalf;
04442 if (ystst < fIMPRy[jl-1]) goto L67;
04443 mnrazz(ystar, fPstar, fIMPRy, jh, jl);
04444 goto L50;
04445 L67:
04446 mnrazz(ystst, fPstst, fIMPRy, jh, jl);
04447 goto L50;
04448
04449 L70:
04450 if (ystar >= fIMPRy[jh-1]) goto L73;
04451 jhold = jh;
04452 mnrazz(ystar, fPstar, fIMPRy, jh, jl);
04453 if (jhold != jh) goto L50;
04454
04455 L73:
04456 for (i = 1; i <= fNpar; ++i) {
04457 fPstst[i-1] = fP[i + jh*fMaxpar - fMaxpar-1]*.5 + fPbar[i-1]*.5;
04458 }
04459 mncalf(fPstst, ycalf);
04460 ystst = ycalf;
04461 if (ystst > fIMPRy[jh-1]) goto L30;
04462
04463 if (ystst < fAmin) goto L67;
04464 mnrazz(ystst, fPstst, fIMPRy, jh, jl);
04465 goto L50;
04466
04467 L95:
04468 if (fISW[4] >= 0) {
04469 Printf(" AN IMPROVEMENT ON THE PREVIOUS MINIMUM HAS BEEN FOUND");
04470 }
04471 reg = .1;
04472
04473 L100:
04474 mninex(fX);
04475 Eval(nparx, fGin, fAmin, fU, 4); ++fNfcn;
04476 for (i = 1; i <= fNpar; ++i) {
04477 fDirin[i-1] = reg*fIMPRdsav[i-1];
04478 if (TMath::Abs(fX[i-1] - fXt[i-1]) > fDirin[i-1]) goto L150;
04479 }
04480 goto L230;
04481 L150:
04482 fNfcnmx = fNfcnmx + npfn - fNfcn;
04483 npfn = fNfcn;
04484 mnsimp();
04485 if (fAmin >= fApsi) goto L325;
04486 for (i = 1; i <= fNpar; ++i) {
04487 fDirin[i-1] = fIMPRdsav[i-1]*.1;
04488 if (TMath::Abs(fX[i-1] - fXt[i-1]) > fDirin[i-1]) goto L250;
04489 }
04490 L230:
04491 if (fAmin < fApsi) goto L350;
04492 goto L325;
04493
04494 L250:
04495 fLnewmn = kTRUE;
04496 if (fISW[1] >= 1) {
04497 fISW[1] = 1;
04498 fDcovar = TMath::Max(fDcovar,.5);
04499 } else fDcovar = 1;
04500 fItaur = 0;
04501 fNfcnmx = fNfcnmx + npfn - fNfcn;
04502 fCstatu = "NEW MINIMU";
04503 if (fISW[4] >= 0) {
04504 Printf(" IMPROVE HAS FOUND A TRULY NEW MINIMUM");
04505 Printf(" *************************************");
04506 }
04507 return;
04508
04509 L280:
04510 if (fISW[4] > 0) {
04511 Printf(" COVARIANCE MATRIX WAS NOT POSITIVE-DEFINITE");
04512 }
04513 goto L325;
04514 L300:
04515 fISW[0] = 1;
04516 L325:
04517 for (i = 1; i <= fNpar; ++i) {
04518 fDirin[i-1] = fIMPRdsav[i-1]*.01;
04519 fX[i-1] = fXt[i-1];
04520 }
04521 fAmin = fApsi;
04522 fEDM = sigsav;
04523 L350:
04524 mninex(fX);
04525 if (fISW[4] > 0) {
04526 Printf(" IMPROVE HAS RETURNED TO REGION OF ORIGINAL MINIMUM");
04527 }
04528 fCstatu = "UNCHANGED ";
04529 mnrset(0);
04530 if (fISW[1] < 2) goto L380;
04531 if (loop < nloop && fISW[0] < 1) goto L20;
04532 L380:
04533 if (iswtr >= 0) mnprin(5, fAmin);
04534 fItaur = 0;
04535 }
04536
04537
04538 void TMinuit::mninex(Double_t *pint)
04539 {
04540
04541
04542
04543
04544
04545
04546 Int_t i, j;
04547
04548 for (j = 0; j < fNpar; ++j) {
04549 i = fNexofi[j]-1;
04550 if (fNvarl[i] == 1) {
04551 fU[i] = pint[j];
04552 } else {
04553 fU[i] = fAlim[i] + (TMath::Sin(pint[j]) + 1)*.5*(fBlim[i] - fAlim[i]);
04554 }
04555 }
04556 }
04557
04558
04559 void TMinuit::mninit(Int_t i1, Int_t i2, Int_t i3)
04560 {
04561
04562
04563
04564
04565
04566
04567
04568 Double_t piby2, epsp1, epsbak, epstry, distnn;
04569 Int_t i, idb;
04570
04571
04572 fIsysrd = i1;
04573 fIsyswr = i2;
04574 fIstkwr[0] = fIsyswr;
04575 fNstkwr = 1;
04576 fIsyssa = i3;
04577 fNstkrd = 0;
04578
04579 fCvrsn = "95.03++ ";
04580
04581 fMaxint = fMaxpar;
04582 fMaxext = 2*fMaxpar;
04583 fUndefi = -54321;
04584 fBigedm = 123456;
04585 fCundef = ")UNDEFINED";
04586 fCovmes[0] = "NO ERROR MATRIX ";
04587 fCovmes[1] = "ERR MATRIX APPROXIMATE";
04588 fCovmes[2] = "ERR MATRIX NOT POS-DEF";
04589 fCovmes[3] = "ERROR MATRIX ACCURATE ";
04590
04591 fNblock = 0;
04592 fIcomnd = 0;
04593 fCtitl = fCundef;
04594 fCfrom = "INPUT ";
04595 fNfcn = 0;
04596 fNfcnfr = fNfcn;
04597 fCstatu = "INITIALIZE";
04598 fISW[2] = 0;
04599 fISW[3] = 0;
04600 fISW[4] = 1;
04601
04602
04603
04604 fISW[5] = 0;
04605
04606
04607 for (idb = 0; idb <= 10; ++idb) { fIdbg[idb] = 0; }
04608 fLrepor = kFALSE;
04609 fLwarn = kTRUE;
04610 fLimset = kFALSE;
04611 fLnewmn = kFALSE;
04612 fIstrat = 1;
04613 fItaur = 0;
04614
04615 fNpagwd = 120;
04616 fNpagln = 56;
04617 fNewpag = 1;
04618 if (fISW[5] > 0) {
04619 fNpagwd = 80;
04620 fNpagln = 30;
04621 fNewpag = 0;
04622 }
04623 fUp = 1;
04624 fUpdflt = fUp;
04625
04626 epstry = .5;
04627 for (i = 1; i <= 100; ++i) {
04628 epstry *= .5;
04629 epsp1 = epstry + 1;
04630 mntiny(epsp1, epsbak);
04631 if (epsbak < epstry) goto L35;
04632 }
04633 epstry = 1e-7;
04634 fEpsmac = epstry*4;
04635 Printf(" MNINIT UNABLE TO DETERMINE ARITHMETIC PRECISION. WILL ASSUME:%g",fEpsmac);
04636 L35:
04637 fEpsmac = epstry*8;
04638 fEpsma2 = TMath::Sqrt(fEpsmac)*2;
04639
04640
04641 piby2 = TMath::ATan(1)*2;
04642 distnn = TMath::Sqrt(fEpsma2)*8;
04643 fVlimhi = piby2 - distnn;
04644 fVlimlo = -piby2 + distnn;
04645 mncler();
04646
04647 }
04648
04649
04650 void TMinuit::mnlims()
04651 {
04652
04653
04654
04655
04656
04657
04658 Double_t dxdi, snew;
04659 Int_t kint, i2, newcod, ifx=0, inu;
04660
04661 fCfrom = "SET LIM ";
04662 fNfcnfr = fNfcn;
04663 fCstatu = "NO CHANGE ";
04664 i2 = Int_t(fWord7[0]);
04665 if (i2 > fMaxext || i2 < 0) goto L900;
04666 if (i2 > 0) goto L30;
04667
04668 newcod = 4;
04669 if (fWord7[1] == fWord7[2]) newcod = 1;
04670 for (inu = 1; inu <= fNu; ++inu) {
04671 if (fNvarl[inu-1] <= 0) continue;
04672 if (fNvarl[inu-1] == 1 && newcod == 1) continue;
04673 kint = fNiofex[inu-1];
04674
04675 if (kint <= 0) {
04676 if (fISW[4] >= 0) {
04677 Printf(" LIMITS NOT CHANGED FOR FIXED PARAMETER:%4d",inu);
04678 }
04679 continue;
04680 }
04681 if (newcod == 1) {
04682
04683 if (fISW[4] > 0) {
04684 Printf(" LIMITS REMOVED FROM PARAMETER :%3d",inu);
04685 }
04686 fCstatu = "NEW LIMITS";
04687 mndxdi(fX[kint-1], kint-1, dxdi);
04688 snew = fGstep[kint-1]*dxdi;
04689 fGstep[kint-1] = TMath::Abs(snew);
04690 fNvarl[inu-1] = 1;
04691 } else {
04692
04693 fAlim[inu-1] = TMath::Min(fWord7[1],fWord7[2]);
04694 fBlim[inu-1] = TMath::Max(fWord7[1],fWord7[2]);
04695 if (fISW[4] > 0) {
04696 Printf(" PARAMETER %3d LIMITS SET TO %15.5g%15.5g",inu,fAlim[inu-1],fBlim[inu-1]);
04697 }
04698 fNvarl[inu-1] = 4;
04699 fCstatu = "NEW LIMITS";
04700 fGstep[kint-1] = -.1;
04701 }
04702 }
04703 goto L900;
04704
04705 L30:
04706 if (fNvarl[i2-1] <= 0) {
04707 Printf(" PARAMETER %3d IS NOT VARIABLE.", i2);
04708 goto L900;
04709 }
04710 kint = fNiofex[i2-1];
04711
04712 if (kint == 0) {
04713 Printf(" REQUEST TO CHANGE LIMITS ON FIXED PARAMETER:%3d",i2);
04714 for (ifx = 1; ifx <= fNpfix; ++ifx) {
04715 if (i2 == fIpfix[ifx-1]) goto L92;
04716 }
04717 Printf(" MINUIT BUG IN MNLIMS. SEE F. JAMES");
04718 L92:
04719 ;
04720 }
04721 if (fWord7[1] != fWord7[2]) goto L235;
04722
04723 if (fNvarl[i2-1] != 1) {
04724 if (fISW[4] > 0) {
04725 Printf(" LIMITS REMOVED FROM PARAMETER %2d",i2);
04726 }
04727 fCstatu = "NEW LIMITS";
04728 if (kint <= 0) {
04729 fGsteps[ifx-1] = TMath::Abs(fGsteps[ifx-1]);
04730 } else {
04731 mndxdi(fX[kint-1], kint-1, dxdi);
04732 if (TMath::Abs(dxdi) < .01) dxdi = .01;
04733 fGstep[kint-1] = TMath::Abs(fGstep[kint-1]*dxdi);
04734 fGrd[kint-1] *= dxdi;
04735 }
04736 fNvarl[i2-1] = 1;
04737 } else {
04738 Printf(" NO LIMITS SPECIFIED. PARAMETER %3d IS ALREADY UNLIMITED. NO CHANGE.",i2);
04739 }
04740 goto L900;
04741
04742 L235:
04743 fAlim[i2-1] = TMath::Min(fWord7[1],fWord7[2]);
04744 fBlim[i2-1] = TMath::Max(fWord7[1],fWord7[2]);
04745 fNvarl[i2-1] = 4;
04746 if (fISW[4] > 0) {
04747 Printf(" PARAMETER %3d LIMITS SET TO %15.5g%15.5g",i2,fAlim[i2-1],fBlim[i2-1]);
04748 }
04749 fCstatu = "NEW LIMITS";
04750 if (kint <= 0) fGsteps[ifx-1] = -.1;
04751 else fGstep[kint-1] = -.1;
04752
04753 L900:
04754 if (fCstatu != "NO CHANGE ") {
04755 mnexin(fX);
04756 mnrset(1);
04757 }
04758 }
04759
04760
04761 void TMinuit::mnline(Double_t *start, Double_t fstart, Double_t *step, Double_t slope, Double_t toler)
04762 {
04763
04764
04765
04766
04767
04768
04769
04770
04771
04772
04773
04774
04775
04776
04777
04778
04779 Double_t xpq[12], ypq[12], slam, sdev, coeff[3], denom, flast;
04780 Double_t fvals[3], xvals[3], f1, fvmin, xvmin, ratio, f2, f3, fvmax;
04781 Double_t toler8, toler9, overal, undral, slamin, slamax, slopem;
04782 Int_t i, nparx=0, nvmax=0, nxypt, kk, ipt;
04783 Bool_t ldebug;
04784 TString cmess;
04785 char chpq[13];
04786 int l65, l70, l80;
04787
04788
04789 l65 = 0; l70 = 0; l80 = 0;
04790 ldebug = fIdbg[1] >= 1;
04791
04792 overal = 1e3;
04793 undral = -100;
04794
04795 if (ldebug) {
04796 mninex(&start[0]);
04797 Eval(nparx, fGin, f1, fU, 4); ++fNfcn;
04798 if (f1 != fstart) {
04799 Printf(" MNLINE start point not consistent, F values, parameters=");
04800 for (kk = 1; kk <= fNpar; ++kk) {
04801 Printf(" %14.5e",fX[kk-1]);
04802 }
04803 }
04804 }
04805
04806 fvmin = fstart;
04807 xvmin = 0;
04808 nxypt = 1;
04809 chpq[0] = charal[0];
04810 xpq[0] = 0;
04811 ypq[0] = fstart;
04812
04813 slamin = 0;
04814 for (i = 1; i <= fNpar; ++i) {
04815 if (step[i-1] != 0) {
04816 ratio = TMath::Abs(start[i-1] / step[i-1]);
04817 if (slamin == 0) slamin = ratio;
04818 if (ratio < slamin) slamin = ratio;
04819 }
04820 fX[i-1] = start[i-1] + step[i-1];
04821 }
04822 if (slamin == 0) slamin = fEpsmac;
04823 slamin *= fEpsma2;
04824 nparx = fNpar;
04825
04826 mninex(fX);
04827 Eval(nparx, fGin, f1, fU, 4); ++fNfcn;
04828 ++nxypt;
04829 chpq[nxypt-1] = charal[nxypt-1];
04830 xpq[nxypt-1] = 1;
04831 ypq[nxypt-1] = f1;
04832 if (f1 < fstart) {
04833 fvmin = f1;
04834 xvmin = 1;
04835 }
04836
04837 slam = 1;
04838 toler8 = toler;
04839 slamax = 5;
04840 flast = f1;
04841
04842
04843 do {
04844 denom = (flast - fstart - slope*slam)*2 / (slam*slam);
04845 slam = 1;
04846 if (denom != 0) slam = -slope / denom;
04847 if (slam < 0) slam = slamax;
04848 if (slam > slamax) slam = slamax;
04849 if (slam < toler8) slam = toler8;
04850 if (slam < slamin) {
04851 l80 = 1;
04852 break;
04853 }
04854 if (TMath::Abs(slam - 1) < toler8 && f1 < fstart) {
04855 l70 = 1;
04856 break;
04857 }
04858 if (TMath::Abs(slam - 1) < toler8) slam = toler8 + 1;
04859 if (nxypt >= 12) {
04860 l65 = 1;
04861 break;
04862 }
04863 for (i = 1; i <= fNpar; ++i) { fX[i-1] = start[i-1] + slam*step[i-1]; }
04864 mninex(fX);
04865 nparx = fNpar;
04866 Eval(nparx, fGin, f2, fU, 4); ++fNfcn;
04867 ++nxypt;
04868 chpq[nxypt-1] = charal[nxypt-1];
04869 xpq[nxypt-1] = slam;
04870 ypq[nxypt-1] = f2;
04871 if (f2 < fvmin) {
04872 fvmin = f2;
04873 xvmin = slam;
04874 }
04875 if (fstart == fvmin) {
04876 flast = f2;
04877 toler8 = toler*slam;
04878 overal = slam - toler8;
04879 slamax = overal;
04880 }
04881 } while (fstart == fvmin);
04882
04883 if (!l65 && !l70 && !l80) {
04884
04885 xvals[0] = xpq[0];
04886 fvals[0] = ypq[0];
04887 xvals[1] = xpq[nxypt-2];
04888 fvals[1] = ypq[nxypt-2];
04889 xvals[2] = xpq[nxypt-1];
04890 fvals[2] = ypq[nxypt-1];
04891
04892 do {
04893 slamax = TMath::Max(slamax,TMath::Abs(xvmin)*2);
04894 mnpfit(xvals, fvals, 3, coeff, sdev);
04895 if (coeff[2] <= 0) {
04896 slopem = coeff[2]*2*xvmin + coeff[1];
04897 if (slopem <= 0) slam = xvmin + slamax;
04898 else slam = xvmin - slamax;
04899 } else {
04900 slam = -coeff[1] / (coeff[2]*2);
04901 if (slam > xvmin + slamax) slam = xvmin + slamax;
04902 if (slam < xvmin - slamax) slam = xvmin - slamax;
04903 }
04904 if (slam > 0) {
04905 if (slam > overal)
04906 slam = overal;
04907 else if (slam < undral)
04908 slam = undral;
04909 }
04910
04911
04912 do {
04913 toler9 = TMath::Max(toler8,TMath::Abs(toler8*slam));
04914 for (ipt = 1; ipt <= 3; ++ipt) {
04915 if (TMath::Abs(slam - xvals[ipt-1]) < toler9) {
04916 l70 = 1;
04917 break;
04918 }
04919 }
04920 if (l70) break;
04921
04922 if (nxypt >= 12) {
04923 l65 = 1;
04924 break;
04925 }
04926 for (i = 1; i <= fNpar; ++i) { fX[i-1] = start[i-1] + slam*step[i-1]; }
04927 mninex(fX);
04928 Eval(nparx, fGin, f3, fU, 4); ++fNfcn;
04929 ++nxypt;
04930 chpq[nxypt-1] = charal[nxypt-1];
04931 xpq[nxypt-1] = slam;
04932 ypq[nxypt-1] = f3;
04933
04934 fvmax = fvals[0];
04935 nvmax = 1;
04936 if (fvals[1] > fvmax) {
04937 fvmax = fvals[1];
04938 nvmax = 2;
04939 }
04940 if (fvals[2] > fvmax) {
04941 fvmax = fvals[2];
04942 nvmax = 3;
04943 }
04944
04945 if (f3 >= fvmax) {
04946 if (nxypt >= 12) {
04947 l65 = 1;
04948 break;
04949 }
04950 if (slam > xvmin) overal = TMath::Min(overal,slam - toler8);
04951 if (slam < xvmin) undral = TMath::Max(undral,slam + toler8);
04952 slam = (slam + xvmin)*.5;
04953 }
04954 } while (f3 >= fvmax);
04955
04956
04957 if (l65 || l70) break;
04958
04959 xvals[nvmax-1] = slam;
04960 fvals[nvmax-1] = f3;
04961 if (f3 < fvmin) {
04962 fvmin = f3;
04963 xvmin = slam;
04964 } else {
04965 if (slam > xvmin) overal = TMath::Min(overal,slam - toler8);
04966 if (slam < xvmin) undral = TMath::Max(undral,slam + toler8);
04967 }
04968 } while (nxypt < 12);
04969 }
04970
04971
04972
04973 if (!l70 && !l80) {
04974 cmess = " LINE SEARCH HAS EXHAUSTED THE LIMIT OF FUNCTION CALLS ";
04975 if (ldebug) {
04976 Printf(" MNLINE DEBUG: steps=");
04977 for (kk = 1; kk <= fNpar; ++kk) {
04978 Printf(" %12.4g",step[kk-1]);
04979 }
04980 }
04981 }
04982
04983 if (l70) cmess = " LINE SEARCH HAS ATTAINED TOLERANCE ";
04984 if (l80) cmess = " STEP SIZE AT ARITHMETICALLY ALLOWED MINIMUM";
04985
04986 fAmin = fvmin;
04987 for (i = 1; i <= fNpar; ++i) {
04988 fDirin[i-1] = step[i-1]*xvmin;
04989 fX[i-1] = start[i-1] + fDirin[i-1];
04990 }
04991 mninex(fX);
04992 if (xvmin < 0) {
04993 mnwarn("D", "MNLINE", " LINE MINIMUM IN BACKWARDS DIRECTION");
04994 }
04995 if (fvmin == fstart) {
04996 mnwarn("D", "MNLINE", " LINE SEARCH FINDS NO IMPROVEMENT ");
04997 }
04998 if (ldebug) {
04999 Printf(" AFTER %3d POINTS,%s",nxypt,(const char*)cmess);
05000 mnplot(xpq, ypq, chpq, nxypt, fNpagwd, fNpagln);
05001 }
05002 }
05003
05004
05005 void TMinuit::mnmatu(Int_t kode)
05006 {
05007
05008
05009
05010
05011
05012
05013
05014 Int_t ndex, i, j, m, n, ncoef, nparm, id, it, ix;
05015 Int_t nsofar, ndi, ndj, iso, isw2, isw5;
05016 TString ctemp;
05017
05018 isw2 = fISW[1];
05019 if (isw2 < 1) {
05020 Printf("%s",(const char*)fCovmes[isw2]);
05021 return;
05022 }
05023 if (fNpar == 0) {
05024 Printf(" MNMATU: NPAR=0");
05025 return;
05026 }
05027
05028 if (kode == 1) {
05029 isw5 = fISW[4];
05030 fISW[4] = 2;
05031 mnemat(fP, fMaxint);
05032 if (isw2 < 3) {
05033 Printf("%s",(const char*)fCovmes[isw2]);
05034 }
05035 fISW[4] = isw5;
05036 }
05037
05038 if (fNpar <= 1) return;
05039 mnwerr();
05040
05041 ncoef = (fNpagwd - 19) / 6;
05042 ncoef = TMath::Min(ncoef,20);
05043 nparm = TMath::Min(fNpar,ncoef);
05044 Printf(" PARAMETER CORRELATION COEFFICIENTS ");
05045 ctemp = " NO. GLOBAL";
05046 for (id = 1; id <= nparm; ++id) {
05047 ctemp += Form(" %6d",fNexofi[id-1]);
05048 }
05049 Printf("%s",(const char*)ctemp);
05050 for (i = 1; i <= fNpar; ++i) {
05051 ix = fNexofi[i-1];
05052 ndi = i*(i + 1) / 2;
05053 for (j = 1; j <= fNpar; ++j) {
05054 m = TMath::Max(i,j);
05055 n = TMath::Min(i,j);
05056 ndex = m*(m-1) / 2 + n;
05057 ndj = j*(j + 1) / 2;
05058 fMATUvline[j-1] = fVhmat[ndex-1] / TMath::Sqrt(TMath::Abs(fVhmat[ndi-1]*fVhmat[ndj-1]));
05059 }
05060 nparm = TMath::Min(fNpar,ncoef);
05061 ctemp = Form(" %3d %7.5f ",ix,fGlobcc[i-1]);
05062 for (it = 1; it <= nparm; ++it) {
05063 ctemp += Form(" %6.3f",fMATUvline[it-1]);
05064 }
05065 Printf("%s",(const char*)ctemp);
05066 if (i <= nparm) continue;
05067 ctemp = " ";
05068 for (iso = 1; iso <= 10; ++iso) {
05069 nsofar = nparm;
05070 nparm = TMath::Min(fNpar,nsofar + ncoef);
05071 for (it = nsofar + 1; it <= nparm; ++it) {
05072 ctemp = ctemp + Form(" %6.3f",fMATUvline[it-1]);
05073 }
05074 Printf("%s",(const char*)ctemp);
05075 if (i <= nparm) break;
05076 }
05077 }
05078 if (isw2 < 3) {
05079 Printf(" %s",(const char*)fCovmes[isw2]);
05080 }
05081 }
05082
05083
05084 void TMinuit::mnmigr()
05085 {
05086
05087
05088
05089
05090
05091
05092
05093
05094 Double_t gdel, gami, vlen, dsum, gssq, vsum, d;
05095 Double_t fzero, fs, ri, delgam, rhotol;
05096 Double_t gdgssq, gvg, vgi;
05097 Int_t npfn, ndex, iext, i, j, m, n, npsdf, nparx;
05098 Int_t iswtr, lined2, kk, nfcnmg, nrstrt,iter;
05099 Bool_t ldebug;
05100 Double_t toler = 0.05;
05101
05102 if (fNpar <= 0) return;
05103 if (fAmin == fUndefi) mnamin();
05104 ldebug = kFALSE; if ( fIdbg[4] >= 1) ldebug = kTRUE;
05105 fCfrom = "MIGRAD ";
05106 fNfcnfr = fNfcn;
05107 nfcnmg = fNfcn;
05108 fCstatu = "INITIATE ";
05109 iswtr = fISW[4] - 2*fItaur;
05110 npfn = fNfcn;
05111 nparx = fNpar;
05112 vlen = (Double_t) (fNpar*(fNpar + 1) / 2);
05113 nrstrt = 0;
05114 npsdf = 0;
05115 lined2 = 0;
05116 fISW[3] = -1;
05117 rhotol = fApsi*.001;
05118 if (iswtr >= 1) {
05119 Printf(" START MIGRAD MINIMIZATION. STRATEGY %2d. CONVERGENCE WHEN EDM .LT.%9.2e",fIstrat,rhotol);
05120 }
05121
05122 if (fIstrat < 2 || fISW[1] >= 3) goto L2;
05123
05124 L1:
05125 if (nrstrt > fIstrat) {
05126 fCstatu = "FAILED ";
05127 fISW[3] = -1;
05128 goto L230;
05129 }
05130
05131 mnhess();
05132 mnwerr();
05133 npsdf = 0;
05134 if (fISW[1] >= 1) goto L10;
05135
05136 L2:
05137 mninex(fX);
05138 if (fISW[2] == 1) {
05139 Eval(nparx, fGin, fzero, fU, 2); ++fNfcn;
05140 }
05141 mnderi();
05142 if (fISW[1] >= 1) goto L10;
05143
05144 for (i = 1; i <= fNpar; ++i) {
05145 fMIGRxxs[i-1] = fX[i-1];
05146 fMIGRstep[i-1] = 0;
05147 }
05148
05149 ++lined2;
05150 if (lined2 < (fIstrat + 1)*fNpar) {
05151 for (i = 1; i <= fNpar; ++i) {
05152 if (fG2[i-1] > 0) continue;
05153 if (fGrd[i-1] > 0) fMIGRstep[i-1] = -TMath::Abs(fGstep[i-1]);
05154 else fMIGRstep[i-1] = TMath::Abs(fGstep[i-1]);
05155 gdel = fMIGRstep[i-1]*fGrd[i-1];
05156 fs = fAmin;
05157 mnline(fMIGRxxs, fs, fMIGRstep, gdel, toler);
05158 mnwarn("D", "MNMIGR", "Negative G2 line search");
05159 iext = fNexofi[i-1];
05160 if (ldebug) {
05161 Printf(" Negative G2 line search, param %3d %13.3g%13.3g",iext,fs,fAmin);
05162 }
05163 goto L2;
05164 }
05165 }
05166
05167 for (i = 1; i <= fNpar; ++i) {
05168 ndex = i*(i-1) / 2;
05169 for (j = 1; j <= i-1; ++j) {
05170 ++ndex;
05171 fVhmat[ndex-1] = 0;
05172 }
05173 ++ndex;
05174 if (fG2[i-1] <= 0) fG2[i-1] = 1;
05175 fVhmat[ndex-1] = 2 / fG2[i-1];
05176 }
05177 fDcovar = 1;
05178 if (ldebug) {
05179 Printf(" DEBUG MNMIGR, STARTING MATRIX DIAGONAL, VHMAT=");
05180 for (kk = 1; kk <= Int_t(vlen); ++kk) {
05181 Printf(" %10.2g",fVhmat[kk-1]);
05182 }
05183 }
05184
05185 L10:
05186 ++nrstrt;
05187 if (nrstrt > fIstrat + 1) {
05188 fCstatu = "FAILED ";
05189 goto L230;
05190 }
05191 fs = fAmin;
05192
05193 fEDM = 0;
05194 for (i = 1; i <= fNpar; ++i) {
05195 fMIGRgs[i-1] = fGrd[i-1];
05196 fMIGRxxs[i-1] = fX[i-1];
05197 ndex = i*(i-1) / 2;
05198 for (j = 1; j <= i-1; ++j) {
05199 ++ndex;
05200 fEDM += fMIGRgs[i-1]*fVhmat[ndex-1]*fMIGRgs[j-1];
05201 }
05202 ++ndex;
05203 fEDM += fMIGRgs[i-1]*fMIGRgs[i-1]*.5*fVhmat[ndex-1];
05204 }
05205 fEDM = fEDM*.5*(fDcovar*3 + 1);
05206 if (fEDM < 0) {
05207 mnwarn("W", "MIGRAD", "STARTING MATRIX NOT POS-DEFINITE.");
05208 fISW[1] = 0;
05209 fDcovar = 1;
05210 goto L2;
05211 }
05212 if (fISW[1] == 0) fEDM = fBigedm;
05213 iter = 0;
05214 mninex(fX);
05215 mnwerr();
05216 if (iswtr >= 1) mnprin(3, fAmin);
05217 if (iswtr >= 2) mnmatu(0);
05218
05219 L24:
05220 if (fNfcn - npfn >= fNfcnmx) goto L190;
05221 gdel = 0;
05222 gssq = 0;
05223 for (i = 1; i <= fNpar; ++i) {
05224 ri = 0;
05225 gssq += fMIGRgs[i-1]*fMIGRgs[i-1];
05226 for (j = 1; j <= fNpar; ++j) {
05227 m = TMath::Max(i,j);
05228 n = TMath::Min(i,j);
05229 ndex = m*(m-1) / 2 + n;
05230 ri += fVhmat[ndex-1]*fMIGRgs[j-1];
05231 }
05232 fMIGRstep[i-1] = ri*-.5;
05233 gdel += fMIGRstep[i-1]*fMIGRgs[i-1];
05234 }
05235 if (gssq == 0) {
05236 mnwarn("D", "MIGRAD", " FIRST DERIVATIVES OF FCN ARE ALL ZERO");
05237 goto L300;
05238 }
05239
05240 if (gdel >= 0) {
05241 mnwarn("D", "MIGRAD", " NEWTON STEP NOT DESCENT.");
05242 if (npsdf == 1) goto L1;
05243 mnpsdf();
05244 npsdf = 1;
05245 goto L24;
05246 }
05247
05248 mnline(fMIGRxxs, fs, fMIGRstep, gdel, toler);
05249 if (fAmin == fs) goto L200;
05250 fCfrom = "MIGRAD ";
05251 fNfcnfr = nfcnmg;
05252 fCstatu = "PROGRESS ";
05253
05254 mninex(fX);
05255 if (fISW[2] == 1) {
05256 Eval(nparx, fGin, fzero, fU, 2); ++fNfcn;
05257 }
05258 mnderi();
05259
05260 npsdf = 0;
05261 L81:
05262 fEDM = 0;
05263 gvg = 0;
05264 delgam = 0;
05265 gdgssq = 0;
05266 for (i = 1; i <= fNpar; ++i) {
05267 ri = 0;
05268 vgi = 0;
05269 for (j = 1; j <= fNpar; ++j) {
05270 m = TMath::Max(i,j);
05271 n = TMath::Min(i,j);
05272 ndex = m*(m-1) / 2 + n;
05273 vgi += fVhmat[ndex-1]*(fGrd[j-1] - fMIGRgs[j-1]);
05274 ri += fVhmat[ndex-1]*fGrd[j-1];
05275 }
05276 fMIGRvg[i-1] = vgi*.5;
05277 gami = fGrd[i-1] - fMIGRgs[i-1];
05278 gdgssq += gami*gami;
05279 gvg += gami*fMIGRvg[i-1];
05280 delgam += fDirin[i-1]*gami;
05281 fEDM += fGrd[i-1]*ri*.5;
05282 }
05283 fEDM = fEDM*.5*(fDcovar*3 + 1);
05284
05285 if (fEDM < 0 || gvg <= 0) {
05286 mnwarn("D", "MIGRAD", "NOT POS-DEF. EDM OR GVG NEGATIVE.");
05287 fCstatu = "NOT POSDEF";
05288 if (npsdf == 1) goto L230;
05289 mnpsdf();
05290 npsdf = 1;
05291 goto L81;
05292 }
05293
05294 ++iter;
05295 if (iswtr >= 3 || (iswtr == 2 && iter % 10 == 1)) {
05296 mnwerr();
05297 mnprin(3, fAmin);
05298 }
05299 if (gdgssq == 0) {
05300 mnwarn("D", "MIGRAD", "NO CHANGE IN FIRST DERIVATIVES OVER LAST STEP");
05301 }
05302 if (delgam < 0) {
05303 mnwarn("D", "MIGRAD", "FIRST DERIVATIVES INCREASING ALONG SEARCH LINE");
05304 }
05305
05306 fCstatu = "IMPROVEMNT";
05307 if (ldebug) {
05308 Printf(" VHMAT 1 =");
05309 for (kk = 1; kk <= 10; ++kk) {
05310 Printf(" %10.2g",fVhmat[kk-1]);
05311 }
05312 }
05313 dsum = 0;
05314 vsum = 0;
05315 for (i = 1; i <= fNpar; ++i) {
05316 for (j = 1; j <= i; ++j) {
05317 if(delgam == 0 || gvg == 0) d = 0;
05318 else d = fDirin[i-1]*fDirin[j-1] / delgam - fMIGRvg[i-1]*fMIGRvg[j-1] / gvg;
05319 dsum += TMath::Abs(d);
05320 ndex = i*(i-1) / 2 + j;
05321 fVhmat[ndex-1] += d*2;
05322 vsum += TMath::Abs(fVhmat[ndex-1]);
05323 }
05324 }
05325
05326 fDcovar = (fDcovar + dsum / vsum)*.5;
05327 if (iswtr >= 3 || ldebug) {
05328 Printf(" RELATIVE CHANGE IN COV. MATRIX=%5.1f per cent",fDcovar*100);
05329 }
05330 if (ldebug) {
05331 Printf(" VHMAT 2 =");
05332 for (kk = 1; kk <= 10; ++kk) {
05333 Printf(" %10.3g",fVhmat[kk-1]);
05334 }
05335 }
05336 if (delgam <= gvg) goto L135;
05337 for (i = 1; i <= fNpar; ++i) {
05338 fMIGRflnu[i-1] = fDirin[i-1] / delgam - fMIGRvg[i-1] / gvg;
05339 }
05340 for (i = 1; i <= fNpar; ++i) {
05341 for (j = 1; j <= i; ++j) {
05342 ndex = i*(i-1) / 2 + j;
05343 fVhmat[ndex-1] += gvg*2*fMIGRflnu[i-1]*fMIGRflnu[j-1];
05344 }
05345 }
05346 L135:
05347
05348 if (fEDM < rhotol*.1) goto L300;
05349
05350 for (i = 1; i <= fNpar; ++i) {
05351 fMIGRxxs[i-1] = fX[i-1];
05352 fMIGRgs[i-1] = fGrd[i-1];
05353 }
05354 fs = fAmin;
05355 if (fISW[1] == 0 && fDcovar < .5) fISW[1] = 1;
05356 if (fISW[1] == 3 && fDcovar > .1) fISW[1] = 1;
05357 if (fISW[1] == 1 && fDcovar < .05) fISW[1] = 3;
05358 goto L24;
05359
05360
05361 L190:
05362 fISW[0] = 1;
05363 if (fISW[4] >= 0) {
05364 Printf(" CALL LIMIT EXCEEDED IN MIGRAD.");
05365 }
05366 fCstatu = "CALL LIMIT";
05367 goto L230;
05368
05369 L200:
05370 if (iswtr >= 1) {
05371 Printf(" MIGRAD FAILS TO FIND IMPROVEMENT");
05372 }
05373 for (i = 1; i <= fNpar; ++i) { fX[i-1] = fMIGRxxs[i-1]; }
05374 if (fEDM < rhotol) goto L300;
05375 if (fEDM < TMath::Abs(fEpsma2*fAmin)) {
05376 if (iswtr >= 0) {
05377 Printf(" MACHINE ACCURACY LIMITS FURTHER IMPROVEMENT.");
05378 }
05379 goto L300;
05380 }
05381 if (fIstrat < 1) {
05382 if (fISW[4] >= 0) {
05383 Printf(" MIGRAD FAILS WITH STRATEGY=0. WILL TRY WITH STRATEGY=1.");
05384 }
05385 fIstrat = 1;
05386 }
05387 goto L1;
05388
05389 L230:
05390 if (iswtr >= 0) {
05391 Printf(" MIGRAD TERMINATED WITHOUT CONVERGENCE.");
05392 }
05393 if (fISW[1] == 3) fISW[1] = 1;
05394 fISW[3] = -1;
05395 goto L400;
05396
05397 L300:
05398 if (iswtr >= 0) {
05399 Printf(" MIGRAD MINIMIZATION HAS CONVERGED.");
05400 }
05401 if (fItaur == 0) {
05402 if (fIstrat >= 2 || (fIstrat == 1 && fISW[1] < 3)) {
05403 if (fISW[4] >= 0) {
05404 Printf(" MIGRAD WILL VERIFY CONVERGENCE AND ERROR MATRIX.");
05405 }
05406 mnhess();
05407 mnwerr();
05408 npsdf = 0;
05409 if (fEDM > rhotol) goto L10;
05410 }
05411 }
05412 fCstatu = "CONVERGED ";
05413 fISW[3] = 1;
05414
05415 L400:
05416 fCfrom = "MIGRAD ";
05417 fNfcnfr = nfcnmg;
05418 mninex(fX);
05419 mnwerr();
05420 if (iswtr >= 0) mnprin(3, fAmin);
05421 if (iswtr >= 1) mnmatu(1);
05422 }
05423
05424
05425 void TMinuit::mnmnos()
05426 {
05427
05428
05429
05430
05431
05432
05433
05434
05435 Double_t val2mi, val2pl;
05436 Int_t nbad, ilax, ilax2, ngood, nfcnmi, iin, knt;
05437
05438 if (fNpar <= 0) goto L700;
05439 ngood = 0;
05440 nbad = 0;
05441 nfcnmi = fNfcn;
05442
05443 for (knt = 1; knt <= fNpar; ++knt) {
05444 if (Int_t(fWord7[1]) == 0) {
05445 ilax = fNexofi[knt-1];
05446 } else {
05447 if (knt >= 7) break;
05448 ilax = Int_t(fWord7[knt]);
05449 if (ilax == 0) break;
05450 if (ilax > 0 && ilax <= fNu) {
05451 if (fNiofex[ilax-1] > 0) goto L565;
05452 }
05453 Printf(" PARAMETER NUMBER %3d NOT A VARIABLE. IGNORED.",ilax);
05454 continue;
05455 }
05456 L565:
05457
05458 ilax2 = 0;
05459 mnmnot(ilax, ilax2, val2pl, val2mi);
05460 if (fLnewmn) goto L650;
05461
05462 iin = fNiofex[ilax-1];
05463 if (fErp[iin-1] > 0) ++ngood;
05464 else ++nbad;
05465 if (fErn[iin-1] < 0) ++ngood;
05466 else ++nbad;
05467 }
05468
05469
05470 fCfrom = "MINOS ";
05471 fNfcnfr = nfcnmi;
05472 fCstatu = "UNCHANGED ";
05473 if (ngood == 0 && nbad == 0) goto L700;
05474 if (ngood > 0 && nbad == 0) fCstatu = "SUCCESSFUL";
05475 if (ngood == 0 && nbad > 0) fCstatu = "FAILURE ";
05476 if (ngood > 0 && nbad > 0) fCstatu = "PROBLEMS ";
05477 if (fISW[4] >= 0) mnprin(4, fAmin);
05478 if (fISW[4] >= 2) mnmatu(0);
05479 return;
05480
05481 L650:
05482 fCfrom = "MINOS ";
05483 fNfcnfr = nfcnmi;
05484 fCstatu = "NEW MINIMU";
05485 if (fISW[4] >= 0) mnprin(4, fAmin);
05486 Printf(" NEW MINIMUM FOUND. GO BACK TO MINIMIZATION STEP.");
05487 Printf(" =================================================");
05488 Printf(" V");
05489 Printf(" V");
05490 Printf(" V");
05491 Printf(" VVVVVVV");
05492 Printf(" VVVVV");
05493 Printf(" VVV");
05494 Printf(" V\n");
05495 return;
05496 L700:
05497 Printf(" THERE ARE NO MINOS ERRORS TO CALCULATE.");
05498 }
05499
05500
05501 void TMinuit::mnmnot(Int_t ilax, Int_t ilax2, Double_t &val2pl, Double_t &val2mi)
05502 {
05503
05504
05505
05506
05507
05508
05509
05510
05511 Int_t i__1;
05512
05513
05514 Double_t delu, aopt, eros;
05515 Double_t abest, xunit, dc, ut, sigsav, du1;
05516 Double_t fac, sig, sav;
05517 Int_t marc, isig, mpar, ndex, imax, indx, ierr, i, j;
05518 Int_t iercr, it, istrav, nfmxin, nlimit, isw2, isw4;
05519 TString csig;
05520
05521
05522 isw2 = fISW[1];
05523 isw4 = fISW[3];
05524 sigsav = fEDM;
05525 istrav = fIstrat;
05526 dc = fDcovar;
05527 fLnewmn = kFALSE;
05528 fApsi = fEpsi*.5;
05529 abest = fAmin;
05530 mpar = fNpar;
05531 nfmxin = fNfcnmx;
05532 for (i = 1; i <= mpar; ++i) { fXt[i-1] = fX[i-1]; }
05533 i__1 = mpar*(mpar + 1) / 2;
05534 for (j = 1; j <= i__1; ++j) { fVthmat[j-1] = fVhmat[j-1]; }
05535 for (i = 1; i <= mpar; ++i) {
05536 fMNOTgcc[i-1] = fGlobcc[i-1];
05537 fMNOTw[i-1] = fWerr[i-1];
05538 }
05539 it = fNiofex[ilax-1];
05540 fErp[it-1] = 0;
05541 fErn[it-1] = 0;
05542 mninex(fXt);
05543 ut = fU[ilax-1];
05544 if (fNvarl[ilax-1] == 1) {
05545 fAlim[ilax-1] = ut - fMNOTw[it-1]*100;
05546 fBlim[ilax-1] = ut + fMNOTw[it-1]*100;
05547 }
05548 ndex = it*(it + 1) / 2;
05549 xunit = TMath::Sqrt(fUp / fVthmat[ndex-1]);
05550 marc = 0;
05551 for (i = 1; i <= mpar; ++i) {
05552 if (i == it) continue;
05553 ++marc;
05554 imax = TMath::Max(it,i);
05555 indx = imax*(imax-1) / 2 + TMath::Min(it,i);
05556 fMNOTxdev[marc-1] = xunit*fVthmat[indx-1];
05557 }
05558
05559 mnfixp(it-1, ierr);
05560 if (ierr > 0) {
05561 Printf(" MINUIT ERROR. CANNOT FIX PARAMETER %4d INTERNAL %3d",ilax,it);
05562 goto L700;
05563 }
05564
05565
05566
05567 for (isig = 1; isig <= 2; ++isig) {
05568 if (isig == 1) {
05569 sig = 1;
05570 csig = "POSI";
05571 } else {
05572 sig = -1;
05573 csig = "NEGA";
05574 }
05575
05576 if (fISW[4] > 1) {
05577 Printf(" DETERMINATION OF %sTIVE MINOS ERROR FOR PARAMETER %d %s"
05578 ,(const char*)csig,ilax
05579 ,(const char*)fCpnam[ilax-1]);
05580 }
05581 if (fISW[1] <= 0) {
05582 mnwarn("D", "MINOS", "NO COVARIANCE MATRIX.");
05583 }
05584 nlimit = fNfcn + nfmxin;
05585 fIstrat = TMath::Max(istrav-1,0);
05586 du1 = fMNOTw[it-1];
05587 fU[ilax-1] = ut + sig*du1;
05588 fU[ilax-1] = TMath::Min(fU[ilax-1],fBlim[ilax-1]);
05589 fU[ilax-1] = TMath::Max(fU[ilax-1],fAlim[ilax-1]);
05590 delu = fU[ilax-1] - ut;
05591
05592 if (TMath::Abs(delu) / (TMath::Abs(ut) + TMath::Abs(fU[ilax-1])) < fEpsmac) goto L440;
05593 fac = delu / fMNOTw[it-1];
05594 for (i = 1; i <= fNpar; ++i) {
05595 fX[i-1] = fXt[i-1] + fac*fMNOTxdev[i-1];
05596 }
05597 if (fISW[4] > 1) {
05598 Printf(" PARAMETER %4d SET TO%11.3e + %10.3e = %12.3e",ilax,ut,delu,fU[ilax-1]);
05599 }
05600
05601 fKe1cr = ilax;
05602 fKe2cr = 0;
05603 fXmidcr = fU[ilax-1];
05604 fXdircr = delu;
05605
05606 fAmin = abest;
05607 fNfcnmx = nlimit - fNfcn;
05608 mncros(aopt, iercr);
05609 if (abest - fAmin > fUp*.01) goto L650;
05610 if (iercr == 1) goto L440;
05611 if (iercr == 2) goto L450;
05612 if (iercr == 3) goto L460;
05613
05614 eros = fXmidcr - ut + aopt*fXdircr;
05615 if (fISW[4] > 1) {
05616 Printf(" THE %4sTIVE MINOS ERROR OF PARAMETER %3d %10s, IS %12.4e"
05617 ,(const char*)csig,ilax
05618 ,(const char*)fCpnam[ilax-1],eros);
05619 }
05620 goto L480;
05621
05622 L440:
05623 if (fISW[4] >= 1) {
05624 Printf(" THE %4sTIVE MINOS ERROR OF PARAMETER %3d, %s EXCEEDS ITS LIMIT."
05625 ,(const char*)csig,ilax
05626 ,(const char*)fCpnam[ilax-1]);
05627 }
05628 eros = fUndefi;
05629 goto L480;
05630 L450:
05631 if (fISW[4] >= 1) {
05632 Printf(" THE %4sTIVE MINOS ERROR %4d REQUIRES MORE THAN %5d FUNCTION CALLS."
05633 ,(const char*)csig,ilax,nfmxin);
05634 }
05635 eros = 0;
05636 goto L480;
05637 L460:
05638 if (fISW[4] >= 1) {
05639 Printf(" %4sTIVE MINOS ERROR NOT CALCULATED FOR PARAMETER %d"
05640 ,(const char*)csig,ilax);
05641 }
05642 eros = 0;
05643
05644 L480:
05645 if (fISW[4] > 1) {
05646 Printf(" **************************************************************************");
05647 }
05648 if (sig < 0) {
05649 fErn[it-1] = eros;
05650 if (ilax2 > 0 && ilax2 <= fNu) val2mi = fU[ilax2-1];
05651 } else {
05652 fErp[it-1] = eros;
05653 if (ilax2 > 0 && ilax2 <= fNu) val2pl = fU[ilax2-1];
05654 }
05655 }
05656
05657
05658 fItaur = 1;
05659 mnfree(1);
05660 i__1 = mpar*(mpar + 1) / 2;
05661 for (j = 1; j <= i__1; ++j) { fVhmat[j-1] = fVthmat[j-1]; }
05662 for (i = 1; i <= mpar; ++i) {
05663 fWerr[i-1] = fMNOTw[i-1];
05664 fGlobcc[i-1] = fMNOTgcc[i-1];
05665 fX[i-1] = fXt[i-1];
05666 }
05667 mninex(fX);
05668 fEDM = sigsav;
05669 fAmin = abest;
05670 fISW[1] = isw2;
05671 fISW[3] = isw4;
05672 fDcovar = dc;
05673 goto L700;
05674
05675 L650:
05676 fLnewmn = kTRUE;
05677 fISW[1] = 0;
05678 fDcovar = 1;
05679 fISW[3] = 0;
05680 sav = fU[ilax-1];
05681 fItaur = 1;
05682 mnfree(1);
05683 fU[ilax-1] = sav;
05684 mnexin(fX);
05685 fEDM = fBigedm;
05686
05687 L700:
05688 fItaur = 0;
05689 fNfcnmx = nfmxin;
05690 fIstrat = istrav;
05691 }
05692
05693
05694 void TMinuit::mnparm(Int_t k1, TString cnamj, Double_t uk, Double_t wk, Double_t a, Double_t b, Int_t &ierflg)
05695 {
05696
05697
05698
05699
05700
05701
05702
05703
05704
05705
05706
05707
05708
05709
05710
05711 Double_t vplu, a_small, gsmin, pinti, vminu, danger, sav, sav2;
05712 Int_t ierr, kint, in, ix, ktofix, lastin, kinfix, nvl;
05713 TString cnamk, chbufi;
05714
05715 Int_t k = k1+1;
05716 cnamk = cnamj;
05717 kint = fNpar;
05718 if (k < 1 || k > fMaxext) {
05719
05720 Printf(" MINUIT USER ERROR. PARAMETER NUMBER IS %3d ALLOWED RANGE IS ONE TO %4d",k,fMaxext);
05721 goto L800;
05722 }
05723
05724 ktofix = 0;
05725 if (fNvarl[k-1] < 0) goto L50;
05726
05727
05728 for (ix = 1; ix <= fNpfix; ++ix) {
05729 if (fIpfix[ix-1] == k) ktofix = k;
05730 }
05731 if (ktofix > 0) {
05732 mnwarn("W", "PARAM DEF", "REDEFINING A FIXED PARAMETER.");
05733 if (kint >= fMaxint) {
05734 Printf(" CANNOT RELEASE. MAX NPAR EXCEEDED.");
05735 goto L800;
05736 }
05737 mnfree(-k);
05738 }
05739
05740 if (fNiofex[k-1] > 0) kint = fNpar - 1;
05741 L50:
05742
05743
05744 if (fLphead && fISW[4] >= 0) {
05745 Printf(" PARAMETER DEFINITIONS:");
05746 Printf(" NO. NAME VALUE STEP SIZE LIMITS");
05747 fLphead = kFALSE;
05748 }
05749 if (wk > 0) goto L122;
05750
05751 if (fISW[4] >= 0) {
05752 Printf(" %5d %-10s %13.5e constant",k,(const char*)cnamk,uk);
05753 }
05754 nvl = 0;
05755 goto L200;
05756 L122:
05757 if (a == 0 && b == 0) {
05758
05759 nvl = 1;
05760 if (fISW[4] >= 0) {
05761 Printf(" %5d %-10s %13.5e%13.5e no limits",k,(const char*)cnamk,uk,wk);
05762 }
05763 } else {
05764
05765 nvl = 4;
05766 fLnolim = kFALSE;
05767 if (fISW[4] >= 0) {
05768 Printf(" %5d %-10s %13.5e%13.5e %13.5e%13.5e",k,(const char*)cnamk,uk,wk,a,b);
05769 }
05770 }
05771
05772 ++kint;
05773 if (kint > fMaxint) {
05774 Printf(" MINUIT USER ERROR. TOO MANY VARIABLE PARAMETERS.");
05775 goto L800;
05776 }
05777 if (nvl == 1) goto L200;
05778 if (a == b) {
05779 Printf(" USER ERROR IN MINUIT PARAMETER");
05780 Printf(" DEFINITION");
05781 Printf(" UPPER AND LOWER LIMITS EQUAL.");
05782 goto L800;
05783 }
05784 if (b < a) {
05785 sav = b;
05786 b = a;
05787 a = sav;
05788 mnwarn("W", "PARAM DEF", "PARAMETER LIMITS WERE REVERSED.");
05789 if (fLwarn) fLphead = kTRUE;
05790 }
05791 if (b - a > 1e7) {
05792 mnwarn("W", "PARAM DEF", Form("LIMITS ON PARAM%d TOO FAR APART.",k));
05793 if (fLwarn) fLphead = kTRUE;
05794 }
05795 danger = (b - uk)*(uk - a);
05796 if (danger < 0) {
05797 mnwarn("W", "PARAM DEF", "STARTING VALUE OUTSIDE LIMITS.");
05798 }
05799 if (danger == 0) {
05800 mnwarn("W", "PARAM DEF", "STARTING VALUE IS AT LIMIT.");
05801 }
05802 L200:
05803
05804
05805 fCfrom = "PARAMETR";
05806 fNfcnfr = fNfcn;
05807 fCstatu = "NEW VALUES";
05808 fNu = TMath::Max(fNu,k);
05809 fCpnam[k-1] = cnamk;
05810 fU[k-1] = uk;
05811 fAlim[k-1] = a;
05812 fBlim[k-1] = b;
05813 fNvarl[k-1] = nvl;
05814 mnrset(1);
05815
05816
05817 lastin = 0;
05818 for (ix = 1; ix <= k-1; ++ix) { if (fNiofex[ix-1] > 0) ++lastin; }
05819
05820 if (kint == fNpar) goto L280;
05821 if (kint > fNpar) {
05822
05823 for (in = fNpar; in >= lastin + 1; --in) {
05824 ix = fNexofi[in-1];
05825 fNiofex[ix-1] = in + 1;
05826 fNexofi[in] = ix;
05827 fX[in] = fX[in-1];
05828 fXt[in] = fXt[in-1];
05829 fDirin[in] = fDirin[in-1];
05830 fG2[in] = fG2[in-1];
05831 fGstep[in] = fGstep[in-1];
05832 fWerr[in] = fWerr[in-1];
05833 fGrd[in] = fGrd[in-1];
05834 }
05835 } else {
05836
05837 for (in = lastin + 1; in <= kint; ++in) {
05838 ix = fNexofi[in];
05839 fNiofex[ix-1] = in;
05840 fNexofi[in-1] = ix;
05841 fX[in-1] = fX[in];
05842 fXt[in-1] = fXt[in];
05843 fDirin[in-1] = fDirin[in];
05844 fG2[in-1] = fG2[in];
05845 fGstep[in-1] = fGstep[in];
05846 fWerr[in-1] = fWerr[in];
05847 fGrd[in-1] = fGrd[in];
05848 }
05849 }
05850 L280:
05851 ix = k;
05852 fNiofex[ix-1] = 0;
05853 fNpar = kint;
05854
05855 if (nvl > 0) {
05856 in = lastin + 1;
05857 fNexofi[in-1] = ix;
05858 fNiofex[ix-1] = in;
05859 sav = fU[ix-1];
05860 mnpint(sav, ix-1, pinti);
05861 fX[in-1] = pinti;
05862 fXt[in-1] = fX[in-1];
05863 fWerr[in-1] = wk;
05864 sav2 = sav + wk;
05865 mnpint(sav2, ix-1, pinti);
05866 vplu = pinti - fX[in-1];
05867 sav2 = sav - wk;
05868 mnpint(sav2, ix-1, pinti);
05869 vminu = pinti - fX[in-1];
05870 fDirin[in-1] = (TMath::Abs(vplu) + TMath::Abs(vminu))*.5;
05871 fG2[in-1] = fUp*2 / (fDirin[in-1]*fDirin[in-1]);
05872 gsmin = fEpsma2*8*TMath::Abs(fX[in-1]);
05873 fGstep[in-1] = TMath::Max(gsmin,fDirin[in-1]*.1);
05874 if (fAmin != fUndefi) {
05875 a_small = TMath::Sqrt(fEpsma2*(fAmin + fUp) / fUp);
05876 fGstep[in-1] = TMath::Max(gsmin,a_small*fDirin[in-1]);
05877 }
05878 fGrd[in-1] = fG2[in-1]*fDirin[in-1];
05879
05880 if (fNvarl[k-1] > 1) {
05881 if (fGstep[in-1] > .5) fGstep[in-1] = .5;
05882 fGstep[in-1] = -fGstep[in-1];
05883 }
05884 }
05885 if (ktofix > 0) {
05886 kinfix = fNiofex[ktofix-1];
05887 if (kinfix > 0) mnfixp(kinfix-1, ierr);
05888 if (ierr > 0) goto L800;
05889 }
05890 ierflg = 0;
05891 return;
05892
05893 L800:
05894 ierflg = 1;
05895 }
05896
05897
05898 void TMinuit::mnpars(TString &crdbuf, Int_t &icondn)
05899 {
05900
05901
05902
05903
05904
05905
05906
05907
05908
05909
05910
05911
05912
05913 Double_t a=0, b=0, fk=0, uk=0, wk=0, xk=0;
05914 Int_t ierr, kapo1, kapo2;
05915 Int_t k, llist, ibegin, lenbuf, istart, lnc, icy;
05916 TString cnamk, comand, celmnt, ctemp;
05917 char stmp[128];
05918
05919 lenbuf = strlen((const char*)crdbuf);
05920
05921 kapo1 = strspn((const char*)crdbuf, "'");
05922 if (kapo1 == 0) goto L150;
05923 kapo2 = strspn((const char*)crdbuf + kapo1, "'");
05924 if (kapo2 == 0) goto L150;
05925
05926 kapo2 += kapo1;
05927
05928 for (istart = 1; istart <= kapo1-1; ++istart) {
05929 if (crdbuf(istart-1,1) != ' ') goto L120;
05930 }
05931 goto L210;
05932 L120:
05933
05934 celmnt = crdbuf(istart-1, kapo1-istart);
05935 if (scanf((const char*)celmnt,&fk)) {;}
05936 k = Int_t(fk);
05937 if (k <= 0) goto L210;
05938 cnamk = "PARAM " + celmnt;
05939 if (kapo2 - kapo1 > 1) {
05940 cnamk = crdbuf(kapo1, kapo2-1-kapo1);
05941 }
05942
05943 for (icy = kapo2 + 1; icy <= lenbuf; ++icy) {
05944 if (crdbuf(icy-1,1) == ',') goto L139;
05945 if (crdbuf(icy-1,1) != ' ') goto L140;
05946 }
05947 uk = 0;
05948 wk = 0;
05949 a = 0;
05950 b = 0;
05951 goto L170;
05952 L139:
05953 ++icy;
05954 L140:
05955 ibegin = icy;
05956 ctemp = crdbuf(ibegin-1,lenbuf-ibegin);
05957 mncrck(ctemp, 20, comand, lnc, fMaxpar, fPARSplist, llist, ierr, fIsyswr);
05958 if (ierr > 0) goto L180;
05959 uk = fPARSplist[0];
05960 wk = 0;
05961 if (llist >= 2) wk = fPARSplist[1];
05962 a = 0;
05963 if (llist >= 3) a = fPARSplist[2];
05964 b = 0;
05965 if (llist >= 4) b = fPARSplist[3];
05966 goto L170;
05967
05968 L150:
05969 if (scanf((const char*)crdbuf,&xk,stmp,&uk,&wk,&a,&b)) {;}
05970 cnamk = stmp;
05971 k = Int_t(xk);
05972 if (k == 0) goto L210;
05973
05974 L170:
05975 mnparm(k-1, cnamk, uk, wk, a, b, ierr);
05976 icondn = ierr;
05977 return;
05978
05979 L180:
05980 icondn = 1;
05981 return;
05982
05983 L210:
05984 icondn = 2;
05985 }
05986
05987
05988 void TMinuit::mnpfit(Double_t *parx2p, Double_t *pary2p, Int_t npar2p, Double_t *coef2p, Double_t &sdev2p)
05989 {
05990
05991
05992
05993
05994
05995
05996
05997
05998
05999
06000
06001
06002
06003 Double_t a, f, s, t, y, s2, x2, x3, x4, y2, cz[3], xm, xy, x2y;
06004 x2 = x3 = 0;
06005 Int_t i;
06006
06007
06008 --coef2p;
06009 --pary2p;
06010 --parx2p;
06011
06012
06013 for (i = 1; i <= 3; ++i) { cz[i-1] = 0; }
06014 sdev2p = 0;
06015 if (npar2p < 3) goto L10;
06016 f = (Double_t) (npar2p);
06017
06018 xm = 0;
06019 for (i = 1; i <= npar2p; ++i) { xm += parx2p[i]; }
06020 xm /= f;
06021 x2 = 0;
06022 x3 = 0;
06023 x4 = 0;
06024 y = 0;
06025 y2 = 0;
06026 xy = 0;
06027 x2y = 0;
06028 for (i = 1; i <= npar2p; ++i) {
06029 s = parx2p[i] - xm;
06030 t = pary2p[i];
06031 s2 = s*s;
06032 x2 += s2;
06033 x3 += s*s2;
06034 x4 += s2*s2;
06035 y += t;
06036 y2 += t*t;
06037 xy += s*t;
06038 x2y += s2*t;
06039 }
06040 a = (f*x4 - x2*x2)*x2 - f*(x3*x3);
06041 if (a == 0) goto L10;
06042 cz[2] = (x2*(f*x2y - x2*y) - f*x3*xy) / a;
06043 cz[1] = (xy - x3*cz[2]) / x2;
06044 cz[0] = (y - x2*cz[2]) / f;
06045 if (npar2p == 3) goto L6;
06046 sdev2p = y2 - (cz[0]*y + cz[1]*xy + cz[2]*x2y);
06047 if (sdev2p < 0) sdev2p = 0;
06048 sdev2p /= f - 3;
06049 L6:
06050 cz[0] += xm*(xm*cz[2] - cz[1]);
06051 cz[1] -= xm*2*cz[2];
06052 L10:
06053 for (i = 1; i <= 3; ++i) { coef2p[i] = cz[i-1]; }
06054 }
06055
06056
06057 void TMinuit::mnpint(Double_t &pexti, Int_t i1, Double_t &pinti)
06058 {
06059
06060
06061
06062
06063
06064
06065 Double_t a, alimi, blimi, yy, yy2;
06066 Int_t igo;
06067 TString chbuf2, chbufi;
06068
06069 Int_t i = i1+1;
06070 pinti = pexti;
06071 igo = fNvarl[i-1];
06072 if (igo == 4) {
06073
06074 alimi = fAlim[i-1];
06075 blimi = fBlim[i-1];
06076 yy = (pexti - alimi)*2 / (blimi - alimi) - 1;
06077 yy2 = yy*yy;
06078 if (yy2 >= 1 - fEpsma2) {
06079 if (yy < 0) {
06080 a = fVlimlo;
06081 chbuf2 = " IS AT ITS LOWER ALLOWED LIMIT.";
06082 } else {
06083 a = fVlimhi;
06084 chbuf2 = " IS AT ITS UPPER ALLOWED LIMIT.";
06085 }
06086 pinti = a;
06087 pexti = alimi + (blimi - alimi)*.5*(TMath::Sin(a) + 1);
06088 fLimset = kTRUE;
06089 if (yy2 > 1) chbuf2 = " BROUGHT BACK INSIDE LIMITS.";
06090 mnwarn("W", fCfrom, Form("VARIABLE%d%s",i,chbuf2.Data()));
06091 } else {
06092 pinti = TMath::ASin(yy);
06093 }
06094 }
06095 }
06096
06097
06098 void TMinuit::mnplot(Double_t *xpt, Double_t *ypt, char *chpt, Int_t nxypt, Int_t npagwd, Int_t npagln)
06099 {
06100
06101
06102
06103
06104
06105
06106
06107
06108
06109
06110
06111
06112
06113
06114
06115
06116
06117 if (fGraphicsMode) {
06118 TPluginHandler *h;
06119 if ((h = gROOT->GetPluginManager()->FindHandler("TMinuitGraph"))) {
06120
06121 if (h->LoadPlugin() != -1)
06122 fPlot = (TObject*)h->ExecPlugin(3,nxypt-2,&xpt[2],&ypt[2]);
06123 }
06124 return;
06125 }
06126
06127 static TString cdot = ".";
06128 static TString cslash = "/";
06129
06130
06131 Double_t xmin, ymin, xmax, ymax, savx, savy, yprt;
06132 Double_t bwidx, bwidy, xbest, ybest, ax, ay, bx, by;
06133 Double_t xvalus[12], any, dxx, dyy;
06134 Int_t iten, i, j, k, maxnx, maxny, iquit, ni, linodd;
06135 Int_t nxbest, nybest, km1, ibk, isp1, nx, ny, ks, ix;
06136 TString chmess, ctemp;
06137 Bool_t overpr;
06138 char cline[120];
06139 char chsav, chbest;
06140
06141
06142
06143 maxnx = TMath::Min(npagwd-20,100);
06144 if (maxnx < 10) maxnx = 10;
06145 maxny = npagln;
06146 if (maxny < 10) maxny = 10;
06147 if (nxypt <= 1) return;
06148 xbest = xpt[0];
06149 ybest = ypt[0];
06150 chbest = chpt[0];
06151
06152 km1 = nxypt - 1;
06153 for (i = 1; i <= km1; ++i) {
06154 iquit = 0;
06155 ni = nxypt - i;
06156 for (j = 1; j <= ni; ++j) {
06157 if (ypt[j-1] > ypt[j]) continue;
06158 savx = xpt[j-1];
06159 xpt[j-1] = xpt[j];
06160 xpt[j] = savx;
06161 savy = ypt[j-1];
06162 ypt[j-1] = ypt[j];
06163 ypt[j] = savy;
06164 chsav = chpt[j-1];
06165 chpt[j-1]= chpt[j];
06166 chpt[j] = chsav;
06167 iquit = 1;
06168 }
06169 if (iquit == 0) break;
06170 }
06171
06172 xmax = xpt[0];
06173 xmin = xmax;
06174 for (i = 1; i <= nxypt; ++i) {
06175 if (xpt[i-1] > xmax) xmax = xpt[i-1];
06176 if (xpt[i-1] < xmin) xmin = xpt[i-1];
06177 }
06178 dxx = (xmax - xmin)*.001;
06179 xmax += dxx;
06180 xmin -= dxx;
06181 mnbins(xmin, xmax, maxnx, xmin, xmax, nx, bwidx);
06182 ymax = ypt[0];
06183 ymin = ypt[nxypt-1];
06184 if (ymax == ymin) ymax = ymin + 1;
06185 dyy = (ymax - ymin)*.001;
06186 ymax += dyy;
06187 ymin -= dyy;
06188 mnbins(ymin, ymax, maxny, ymin, ymax, ny, bwidy);
06189 any = (Double_t) ny;
06190
06191 if (chbest == ' ') goto L50;
06192 xbest = (xmax + xmin)*.5;
06193 ybest = (ymax + ymin)*.5;
06194 L50:
06195
06196 ax = 1 / bwidx;
06197 ay = 1 / bwidy;
06198 bx = -ax*xmin + 2;
06199 by = -ay*ymin - 2;
06200
06201 for (i = 1; i <= nxypt; ++i) {
06202 xpt[i-1] = ax*xpt[i-1] + bx;
06203 ypt[i-1] = any - ay*ypt[i-1] - by;
06204 }
06205 nxbest = Int_t((ax*xbest + bx));
06206 nybest = Int_t((any - ay*ybest - by));
06207
06208 ny += 2;
06209 nx += 2;
06210 isp1 = 1;
06211 linodd = 1;
06212 overpr = kFALSE;
06213 for (i = 1; i <= ny; ++i) {
06214 for (ibk = 1; ibk <= nx; ++ibk) { cline[ibk-1] = ' '; }
06215 cline[nx] = '\0';
06216 cline[nx+1] = '\0';
06217 cline[0] = '.';
06218 cline[nx-1] = '.';
06219 cline[nxbest-1] = '.';
06220 if (i != 1 && i != nybest && i != ny) goto L320;
06221 for (j = 1; j <= nx; ++j) { cline[j-1] = '.'; }
06222 L320:
06223 yprt = ymax - Double_t(i-1)*bwidy;
06224 if (isp1 > nxypt) goto L350;
06225
06226 for (k = isp1; k <= nxypt; ++k) {
06227 ks = Int_t(ypt[k-1]);
06228 if (ks > i) goto L345;
06229 ix = Int_t(xpt[k-1]);
06230 if (cline[ix-1] == '.') goto L340;
06231 if (cline[ix-1] == ' ') goto L340;
06232 if (cline[ix-1] == chpt[k-1]) continue;
06233 overpr = kTRUE;
06234
06235
06236 cline[ix-1] = '&';
06237 continue;
06238 L340:
06239 cline[ix-1] = chpt[k-1];
06240 }
06241 isp1 = nxypt + 1;
06242 goto L350;
06243 L345:
06244 isp1 = k;
06245 L350:
06246 if (linodd == 1 || i == ny) goto L380;
06247 linodd = 1;
06248 ctemp = cline;
06249 Printf(" %s",(const char*)ctemp);
06250 goto L400;
06251 L380:
06252 ctemp = cline;
06253 Printf(" %14.7g ..%s",yprt,(const char*)ctemp);
06254 linodd = 0;
06255 L400:
06256 ;
06257 }
06258
06259 for (ibk = 1; ibk <= nx; ++ibk) {
06260 cline[ibk-1] = ' ';
06261 if (ibk % 10 == 1) cline[ibk-1] = '/';
06262 }
06263 Printf(" %s",cline);
06264
06265 for (ibk = 1; ibk <= 12; ++ibk) {
06266 xvalus[ibk-1] = xmin + Double_t(ibk-1)*10*bwidx;
06267 }
06268 printf(" ");
06269 iten = (nx + 9) / 10;
06270 for (ibk = 1; ibk <= iten; ++ibk) {
06271 Printf(" %9.4g", xvalus[ibk-1]);
06272 }
06273 chmess = " ";
06274 if (overpr) chmess = " Overprint character is &";
06275 Printf(" ONE COLUMN=%13.7g%s",bwidx,(const char*)chmess);
06276 }
06277
06278
06279 void TMinuit::mnpout(Int_t iuext1, TString &chnam, Double_t &val, Double_t &err, Double_t &xlolim, Double_t &xuplim, Int_t &iuint) const
06280 {
06281
06282
06283
06284
06285
06286
06287
06288
06289
06290
06291
06292
06293
06294
06295
06296
06297
06298 Int_t iint, iext, nvl;
06299
06300 Int_t iuext = iuext1 + 1;
06301 xlolim = 0;
06302 xuplim = 0;
06303 err = 0;
06304 if (iuext == 0) goto L100;
06305 if (iuext < 0) {
06306
06307 iint = -(iuext);
06308 if (iint > fNpar) goto L100;
06309 iext = fNexofi[iint-1];
06310 iuint = iext;
06311 } else {
06312
06313 iext = iuext;
06314 if (iext > fNu) goto L100;
06315 iint = fNiofex[iext-1];
06316 iuint = iint;
06317 }
06318
06319 nvl = fNvarl[iext-1];
06320 if (nvl < 0) goto L100;
06321 chnam = fCpnam[iext-1];
06322 val = fU[iext-1];
06323 if (iint > 0) err = fWerr[iint-1];
06324 if (nvl == 4) {
06325 xlolim = fAlim[iext-1];
06326 xuplim = fBlim[iext-1];
06327 }
06328 return;
06329
06330 L100:
06331 iuint = -1;
06332 chnam = "undefined";
06333 val = 0;
06334 }
06335
06336
06337 void TMinuit::mnprin(Int_t inkode, Double_t fval)
06338 {
06339
06340
06341
06342
06343
06344
06345
06346
06347
06348
06349
06350
06351
06352
06353
06354
06355 static TString cblank = " ";
06356 static TString cnambf = " ";
06357
06358
06359 Double_t dcmax, x1, x2, x3, dc;
06360 x2 = x3 = 0;
06361 Int_t nadd, i, k, l, m, ikode, ic, nc, ntrail, lbl;
06362 TString chedm;
06363 TString colhdl[6], colhdu[6], cx2, cx3, cheval;
06364
06365 if (fNu == 0) {
06366 Printf(" THERE ARE CURRENTLY NO PARAMETERS DEFINED");
06367 return;
06368 }
06369
06370 ikode = inkode;
06371 if (inkode == 5) {
06372 ikode = fISW[1] + 1;
06373 if (ikode > 3) ikode = 3;
06374 }
06375
06376 for (k = 1; k <= 6; ++k) {
06377 colhdu[k-1] = "UNDEFINED";
06378 colhdl[k-1] = "COLUMN HEAD";
06379 }
06380
06381 if (ikode == 4 && fCtitl != fCundef) {
06382 Printf(" MINUIT TASK: %s",(const char*)fCtitl);
06383 }
06384
06385 if (fval == fUndefi) cheval = " unknown ";
06386 else cheval = Form("%g",fval);
06387
06388 if (fEDM == fBigedm) chedm = " unknown ";
06389 else chedm = Form("%g",fEDM);
06390
06391 nc = fNfcn - fNfcnfr;
06392 Printf(" FCN=%s FROM %8s STATUS=%10s %6d CALLS %9d TOTAL"
06393 ,(const char*)cheval
06394 ,(const char*)fCfrom
06395 ,(const char*)fCstatu,nc,fNfcn);
06396 m = fISW[1];
06397 if (m == 0 || m == 2 || fDcovar == 0) {
06398 Printf(" EDM=%s STRATEGY=%2d %s"
06399 ,(const char*)chedm,fIstrat
06400 ,(const char*)fCovmes[m]);
06401 } else {
06402 dcmax = 1;
06403 dc = TMath::Min(fDcovar,dcmax)*100;
06404 Printf(" EDM=%s STRATEGY=%2d ERROR MATRIX UNCERTAINTY %5.1f per cent"
06405 ,(const char*)chedm,fIstrat,dc);
06406 }
06407
06408 if (ikode == 0) return;
06409
06410 ntrail = 10;
06411 for (i = 1; i <= fNu; ++i) {
06412 if (fNvarl[i-1] < 0) continue;
06413 for (ic = 10; ic >= 1; --ic) {
06414 if (fCpnam[i-1](ic-1,1) != " ") goto L16;
06415 }
06416 ic = 1;
06417 L16:
06418 lbl = 10 - ic;
06419 if (lbl < ntrail) ntrail = lbl;
06420 }
06421 nadd = ntrail / 2 + 1;
06422 if (ikode == 1) {
06423 colhdu[0] = " ";
06424 colhdl[0] = " ERROR ";
06425 colhdu[1] = " PHYSICAL";
06426 colhdu[2] = " LIMITS ";
06427 colhdl[1] = " NEGATIVE ";
06428 colhdl[2] = " POSITIVE ";
06429 }
06430 if (ikode == 2) {
06431 colhdu[0] = " ";
06432 colhdl[0] = " ERROR ";
06433 colhdu[1] = " INTERNAL ";
06434 colhdl[1] = " STEP SIZE ";
06435 colhdu[2] = " INTERNAL ";
06436 colhdl[2] = " VALUE ";
06437 }
06438 if (ikode == 3) {
06439 colhdu[0] = " ";
06440 colhdl[0] = " ERROR ";
06441 colhdu[1] = " STEP ";
06442 colhdl[1] = " SIZE ";
06443 colhdu[2] = " FIRST ";
06444 colhdl[2] = " DERIVATIVE ";
06445 }
06446 if (ikode == 4) {
06447 colhdu[0] = " PARABOLIC ";
06448 colhdl[0] = " ERROR ";
06449 colhdu[1] = " MINOS ";
06450 colhdu[2] = "ERRORS ";
06451 colhdl[1] = " NEGATIVE ";
06452 colhdl[2] = " POSITIVE ";
06453 }
06454
06455 if (ikode != 4) {
06456 if (fISW[1] < 3) colhdu[0] = " APPROXIMATE ";
06457 if (fISW[1] < 1) colhdu[0] = " CURRENT GUESS";
06458 }
06459 Printf(" EXT PARAMETER %-14s%-14s%-14s",(const char*)colhdu[0]
06460 ,(const char*)colhdu[1]
06461 ,(const char*)colhdu[2]);
06462 Printf(" NO. NAME VALUE %-14s%-14s%-14s",(const char*)colhdl[0]
06463 ,(const char*)colhdl[1]
06464 ,(const char*)colhdl[2]);
06465
06466 for (i = 1; i <= fNu; ++i) {
06467 if (fNvarl[i-1] < 0) continue;
06468 l = fNiofex[i-1];
06469 cnambf = cblank(0,nadd) + fCpnam[i-1];
06470 if (l == 0) goto L55;
06471
06472 x1 = fWerr[l-1];
06473 cx2 = "PLEASE GET X..";
06474 cx3 = "PLEASE GET X..";
06475 if (ikode == 1) {
06476 if (fNvarl[i-1] <= 1) {
06477 Printf("%4d %-11s%14.5e%14.5e",i,(const char*)cnambf,fU[i-1],x1);
06478 continue;
06479 } else {
06480 x2 = fAlim[i-1];
06481 x3 = fBlim[i-1];
06482 }
06483 }
06484 if (ikode == 2) {
06485 x2 = fDirin[l-1];
06486 x3 = fX[l-1];
06487 }
06488 if (ikode == 3) {
06489 x2 = fDirin[l-1];
06490 x3 = fGrd[l-1];
06491 if (fNvarl[i-1] > 1 && TMath::Abs(TMath::Cos(fX[l-1])) < .001) {
06492 cx3 = "** at limit **";
06493 }
06494 }
06495 if (ikode == 4) {
06496 x2 = fErn[l-1];
06497 if (x2 == 0) cx2 = " ";
06498 if (x2 == fUndefi) cx2 = " at limit ";
06499 x3 = fErp[l-1];
06500 if (x3 == 0) cx3 = " ";
06501 if (x3 == fUndefi) cx3 = " at limit ";
06502 }
06503 if (cx2 == "PLEASE GET X..") cx2 = Form("%14.5e",x2);
06504 if (cx3 == "PLEASE GET X..") cx3 = Form("%14.5e",x3);
06505 Printf("%4d %-11s%14.5e%14.5e%-14s%-14s",i
06506 ,(const char*)cnambf,fU[i-1],x1
06507 ,(const char*)cx2,(const char*)cx3);
06508
06509
06510 if (fNvarl[i-1] <= 1 || ikode == 3) continue;
06511 if (TMath::Abs(TMath::Cos(fX[l-1])) < .001) {
06512 Printf(" WARNING - - ABOVE PARAMETER IS AT LIMIT.");
06513 }
06514 continue;
06515
06516
06517 L55:
06518 colhdu[0] = " constant ";
06519 if (fNvarl[i-1] > 0) colhdu[0] = " fixed ";
06520 if (fNvarl[i-1] == 4 && ikode == 1) {
06521 Printf("%4d %-11s%14.5e%-14s%14.5e%14.5e",i
06522 ,(const char*)cnambf,fU[i-1]
06523 ,(const char*)colhdu[0],fAlim[i-1],fBlim[i-1]);
06524 } else {
06525 Printf("%4d %-11s%14.5e%s",i
06526 ,(const char*)cnambf,fU[i-1],(const char*)colhdu[0]);
06527 }
06528 }
06529
06530 if (fUp != fUpdflt) {
06531 Printf(" ERR DEF= %g",fUp);
06532 }
06533 return;
06534 }
06535
06536
06537 void TMinuit::mnpsdf()
06538 {
06539
06540
06541
06542
06543
06544
06545 Double_t dgmin, padd, pmin, pmax, dg, epspdf, epsmin;
06546 Int_t ndex, i, j, ndexd, ip, ifault;
06547 TString chbuff, ctemp;
06548
06549 epsmin = 1e-6;
06550 epspdf = TMath::Max(epsmin,fEpsma2);
06551 dgmin = fVhmat[0];
06552
06553 for (i = 1; i <= fNpar; ++i) {
06554 ndex = i*(i + 1) / 2;
06555 if (fVhmat[ndex-1] <= 0) {
06556 mnwarn("W", fCfrom, Form("Negative diagonal element %d in Error Matrix",i));
06557 }
06558 if (fVhmat[ndex-1] < dgmin) dgmin = fVhmat[ndex-1];
06559 }
06560 if (dgmin <= 0) {
06561 dg = epspdf + 1 - dgmin;
06562 mnwarn("W", fCfrom, Form("%g added to diagonal of error matrix",dg));
06563 } else {
06564 dg = 0;
06565 }
06566
06567 for (i = 1; i <= fNpar; ++i) {
06568 ndex = i*(i-1) / 2;
06569 ndexd = ndex + i;
06570 fVhmat[ndexd-1] += dg;
06571 if (fVhmat[ndexd-1]==0) {
06572 fPSDFs[i-1] = 1 / 1e-19;
06573 } else {
06574 fPSDFs[i-1] = 1 / TMath::Sqrt(fVhmat[ndexd-1]);
06575 }
06576 for (j = 1; j <= i; ++j) {
06577 ++ndex;
06578 fP[i + j*fMaxpar - fMaxpar-1] = fVhmat[ndex-1]*fPSDFs[i-1]*fPSDFs[j-1];
06579 }
06580 }
06581
06582 mneig(fP, fMaxint, fNpar, fMaxint, fPstar, epspdf, ifault);
06583 pmin = fPstar[0];
06584 pmax = fPstar[0];
06585 for (ip = 2; ip <= fNpar; ++ip) {
06586 if (fPstar[ip-1] < pmin) pmin = fPstar[ip-1];
06587 if (fPstar[ip-1] > pmax) pmax = fPstar[ip-1];
06588 }
06589 pmax = TMath::Max(TMath::Abs(pmax),Double_t(1));
06590 if ((pmin <= 0 && fLwarn) || fISW[4] >= 2) {
06591 Printf(" EIGENVALUES OF SECOND-DERIVATIVE MATRIX:");
06592 ctemp = " ";
06593 for (ip = 1; ip <= fNpar; ++ip) {
06594 ctemp += Form(" %11.4e",fPstar[ip-1]);
06595 }
06596 Printf("%s", ctemp.Data());
06597 }
06598 if (pmin > epspdf*pmax) return;
06599 if (fISW[1] == 3) fISW[1] = 2;
06600 padd = pmax*.001 - pmin;
06601 for (ip = 1; ip <= fNpar; ++ip) {
06602 ndex = ip*(ip + 1) / 2;
06603 fVhmat[ndex-1] *= padd + 1;
06604 }
06605 fCstatu = "NOT POSDEF";
06606 mnwarn("W", fCfrom, Form("MATRIX FORCED POS-DEF BY ADDING %f TO DIAGONAL.",padd));
06607
06608 }
06609
06610
06611 void TMinuit::mnrazz(Double_t ynew, Double_t *pnew, Double_t *y, Int_t &jh, Int_t &jl)
06612 {
06613
06614
06615
06616
06617
06618
06619
06620 Double_t pbig, plit;
06621 Int_t i, j, nparp1;
06622
06623
06624 for (i = 1; i <= fNpar; ++i) { fP[i + jh*fMaxpar - fMaxpar-1] = pnew[i-1]; }
06625 y[jh-1] = ynew;
06626 if (ynew < fAmin) {
06627 for (i = 1; i <= fNpar; ++i) { fX[i-1] = pnew[i-1]; }
06628 mninex(fX);
06629 fAmin = ynew;
06630 fCstatu = "PROGRESS ";
06631 jl = jh;
06632 }
06633 jh = 1;
06634 nparp1 = fNpar + 1;
06635 for (j = 2; j <= nparp1; ++j) { if (y[j-1] > y[jh-1]) jh = j; }
06636 fEDM = y[jh-1] - y[jl-1];
06637 if (fEDM <= 0) goto L45;
06638 for (i = 1; i <= fNpar; ++i) {
06639 pbig = fP[i-1];
06640 plit = pbig;
06641 for (j = 2; j <= nparp1; ++j) {
06642 if (fP[i + j*fMaxpar - fMaxpar-1] > pbig) pbig = fP[i + j*fMaxpar - fMaxpar-1];
06643 if (fP[i + j*fMaxpar - fMaxpar-1] < plit) plit = fP[i + j*fMaxpar - fMaxpar-1];
06644 }
06645 fDirin[i-1] = pbig - plit;
06646 }
06647 L40:
06648 return;
06649 L45:
06650 Printf(" FUNCTION VALUE DOES NOT SEEM TO DEPEND ON ANY OF THE %d VARIABLE PARAMETERS.",fNpar);
06651 Printf(" VERIFY THAT STEP SIZES ARE BIG ENOUGH AND CHECK FCN LOGIC.");
06652 Printf(" *******************************************************************************");
06653 Printf(" *******************************************************************************");
06654 goto L40;
06655 }
06656
06657
06658 void TMinuit::mnrn15(Double_t &val, Int_t &inseed)
06659 {
06660
06661
06662
06663
06664
06665
06666
06667
06668
06669
06670 static Int_t iseed = 12345;
06671
06672 Int_t k;
06673
06674 if (val == 3) goto L100;
06675 inseed = iseed;
06676 k = iseed / 53668;
06677 iseed = (iseed - k*53668)*40014 - k*12211;
06678 if (iseed < 0) iseed += 2147483563;
06679 val = Double_t(iseed*4.656613e-10);
06680 return;
06681
06682 L100:
06683 iseed = inseed;
06684 }
06685
06686
06687 void TMinuit::mnrset(Int_t iopt)
06688 {
06689
06690
06691
06692
06693
06694
06695
06696
06697 Int_t iext, i;
06698
06699 fCstatu = "RESET ";
06700 if (iopt >= 1) {
06701 fAmin = fUndefi;
06702 fFval3 = TMath::Abs(fAmin)*2 + 1;
06703 fEDM = fBigedm;
06704 fISW[3] = 0;
06705 fISW[1] = 0;
06706 fDcovar = 1;
06707 fISW[0] = 0;
06708 }
06709 fLnolim = kTRUE;
06710 for (i = 1; i <= fNpar; ++i) {
06711 iext = fNexofi[i-1];
06712 if (fNvarl[iext-1] >= 4) fLnolim = kFALSE;
06713 fErp[i-1] = 0;
06714 fErn[i-1] = 0;
06715 fGlobcc[i-1] = 0;
06716 }
06717 if (fISW[1] >= 1) {
06718 fISW[1] = 1;
06719 fDcovar = TMath::Max(fDcovar,.5);
06720 }
06721 }
06722
06723
06724 void TMinuit::mnsave()
06725 {
06726
06727
06728
06729
06730
06731
06732 Printf("mnsave is dummy in TMinuit");
06733
06734 }
06735
06736
06737 void TMinuit::mnscan()
06738 {
06739
06740
06741
06742
06743
06744
06745
06746
06747 Double_t step, uhigh, xhreq, xlreq, ubest, fnext, unext, xh, xl;
06748 Int_t ipar, iint, icall, ncall, nbins, nparx;
06749 Int_t nxypt, nccall, iparwd;
06750
06751 xlreq = TMath::Min(fWord7[2],fWord7[3]);
06752 xhreq = TMath::Max(fWord7[2],fWord7[3]);
06753 ncall = Int_t((fWord7[1] + .01));
06754 if (ncall <= 1) ncall = 41;
06755 if (ncall > 98) ncall = 98;
06756 nccall = ncall;
06757 if (fAmin == fUndefi) mnamin();
06758 iparwd = Int_t((fWord7[0] + .1));
06759 ipar = TMath::Max(iparwd,0);
06760 fCstatu = "NO CHANGE";
06761 if (iparwd > 0) goto L200;
06762
06763
06764 L100:
06765 ++ipar;
06766 if (ipar > fNu) goto L900;
06767 iint = fNiofex[ipar-1];
06768 if (iint <= 0) goto L100;
06769
06770 L200:
06771 iint = fNiofex[ipar-1];
06772 ubest = fU[ipar-1];
06773 fXpt[0] = ubest;
06774 fYpt[0] = fAmin;
06775 fChpt[0] = ' ';
06776 fXpt[1] = ubest;
06777 fYpt[1] = fAmin;
06778 fChpt[1] = 'X';
06779 nxypt = 2;
06780 if (fNvarl[ipar-1] > 1) goto L300;
06781
06782
06783 if (xlreq == xhreq) goto L250;
06784 unext = xlreq;
06785 step = (xhreq - xlreq) / Double_t(ncall-1);
06786 goto L500;
06787 L250:
06788 xl = ubest - fWerr[iint-1];
06789 xh = ubest + fWerr[iint-1];
06790 mnbins(xl, xh, ncall, unext, uhigh, nbins, step);
06791 nccall = nbins + 1;
06792 goto L500;
06793
06794 L300:
06795 if (xlreq == xhreq) goto L350;
06796
06797 xl = TMath::Max(xlreq,fAlim[ipar-1]);
06798
06799 xh = TMath::Min(xhreq,fBlim[ipar-1]);
06800 if (xl >= xh) goto L700;
06801 unext = xl;
06802 step = (xh - xl) / Double_t(ncall-1);
06803 goto L500;
06804 L350:
06805 unext = fAlim[ipar-1];
06806 step = (fBlim[ipar-1] - fAlim[ipar-1]) / Double_t(ncall-1);
06807
06808 L500:
06809 for (icall = 1; icall <= nccall; ++icall) {
06810 fU[ipar-1] = unext;
06811 nparx = fNpar;
06812 Eval(nparx, fGin, fnext, fU, 4); ++fNfcn;
06813 ++nxypt;
06814 fXpt[nxypt-1] = unext;
06815 fYpt[nxypt-1] = fnext;
06816 fChpt[nxypt-1] = '*';
06817 if (fnext < fAmin) {
06818 fAmin = fnext;
06819 ubest = unext;
06820 fCstatu = "IMPROVED ";
06821 }
06822 unext += step;
06823 }
06824 fChpt[nccall] = 0;
06825
06826
06827 fU[ipar-1] = ubest;
06828 mnexin(fX);
06829 if (fISW[4] >= 1)
06830 Printf("%dSCAN OF PARAMETER NO. %d, %s"
06831 ,fNewpag,ipar,(const char*)fCpnam[ipar-1]);
06832 mnplot(fXpt, fYpt, fChpt, nxypt, fNpagwd, fNpagln);
06833 goto L800;
06834 L700:
06835 Printf(" REQUESTED RANGE OUTSIDE LIMITS FOR PARAMETER %d",ipar);
06836 L800:
06837 if (iparwd <= 0) goto L100;
06838
06839 L900:
06840 if (fISW[4] >= 0) mnprin(5, fAmin);
06841 }
06842
06843
06844 void TMinuit::mnseek()
06845 {
06846
06847
06848
06849
06850
06851
06852
06853
06854
06855
06856
06857
06858 Double_t dxdi, rnum, ftry, rnum1, rnum2, alpha;
06859 Double_t flast, bar;
06860 Int_t ipar, iext, j, ifail, iseed=0, nparx, istep, ib, mxfail, mxstep;
06861
06862 mxfail = Int_t(fWord7[0]);
06863 if (mxfail <= 0) mxfail = fNpar*20 + 100;
06864 mxstep = mxfail*10;
06865 if (fAmin == fUndefi) mnamin();
06866 alpha = fWord7[1];
06867 if (alpha <= 0) alpha = 3;
06868 if (fISW[4] >= 1) {
06869 Printf(" MNSEEK: MONTE CARLO MINIMIZATION USING METROPOLIS ALGORITHM");
06870 Printf(" TO STOP AFTER %6d SUCCESSIVE FAILURES, OR %7d STEPS",mxfail,mxstep);
06871 Printf(" MAXIMUM STEP SIZE IS %9.3f ERROR BARS.",alpha);
06872 }
06873 fCstatu = "INITIAL ";
06874 if (fISW[4] >= 2) mnprin(2, fAmin);
06875 fCstatu = "UNCHANGED ";
06876 ifail = 0;
06877 rnum = 0;
06878 rnum1 = 0;
06879 rnum2 = 0;
06880 nparx = fNpar;
06881 flast = fAmin;
06882
06883 for (ipar = 1; ipar <= fNpar; ++ipar) {
06884 iext = fNexofi[ipar-1];
06885 fDirin[ipar-1] = alpha*2*fWerr[ipar-1];
06886 if (fNvarl[iext-1] > 1) {
06887
06888 mndxdi(fX[ipar-1], ipar-1, dxdi);
06889 if (dxdi == 0) dxdi = 1;
06890 fDirin[ipar-1] = alpha*2*fWerr[ipar-1] / dxdi;
06891 if (TMath::Abs(fDirin[ipar-1]) > 6.2831859999999997) {
06892 fDirin[ipar-1] = 6.2831859999999997;
06893 }
06894 }
06895 fSEEKxmid[ipar-1] = fX[ipar-1];
06896 fSEEKxbest[ipar-1] = fX[ipar-1];
06897 }
06898
06899 for (istep = 1; istep <= mxstep; ++istep) {
06900 if (ifail >= mxfail) break;
06901 for (ipar = 1; ipar <= fNpar; ++ipar) {
06902 mnrn15(rnum1, iseed);
06903 mnrn15(rnum2, iseed);
06904 fX[ipar-1] = fSEEKxmid[ipar-1] + (rnum1 + rnum2 - 1)*.5*fDirin[ipar-1];
06905 }
06906 mninex(fX);
06907 Eval(nparx, fGin, ftry, fU, 4); ++fNfcn;
06908 if (ftry < flast) {
06909 if (ftry < fAmin) {
06910 fCstatu = "IMPROVEMNT";
06911 fAmin = ftry;
06912 for (ib = 1; ib <= fNpar; ++ib) { fSEEKxbest[ib-1] = fX[ib-1]; }
06913 ifail = 0;
06914 if (fISW[4] >= 2) mnprin(2, fAmin);
06915 }
06916 goto L300;
06917 } else {
06918 ++ifail;
06919
06920 bar = (fAmin - ftry) / fUp;
06921 mnrn15(rnum, iseed);
06922 if (bar < TMath::Log(rnum)) continue;
06923 }
06924
06925 L300:
06926 for (j = 1; j <= fNpar; ++j) { fSEEKxmid[j-1] = fX[j-1]; }
06927 flast = ftry;
06928 }
06929
06930 if (fISW[4] > 1) {
06931 Printf(" MNSEEK: %5d SUCCESSIVE UNSUCCESSFUL TRIALS.",ifail);
06932 }
06933 for (ib = 1; ib <= fNpar; ++ib) { fX[ib-1] = fSEEKxbest[ib-1]; }
06934 mninex(fX);
06935 if (fISW[4] >= 1) mnprin(2, fAmin);
06936 if (fISW[4] == 0) mnprin(0, fAmin);
06937 }
06938
06939
06940 void TMinuit::mnset()
06941 {
06942
06943
06944
06945
06946
06947
06948
06949
06950
06951
06952
06953
06954
06955 const char *cname[30] = {
06956 "FCN value ",
06957 "PARameters",
06958 "LIMits ",
06959 "COVariance",
06960 "CORrelatio",
06961 "PRInt levl",
06962 "NOGradient",
06963 "GRAdient ",
06964 "ERRor def ",
06965 "INPut file",
06966 "WIDth page",
06967 "LINes page",
06968 "NOWarnings",
06969 "WARnings ",
06970 "RANdom gen",
06971 "TITle ",
06972 "STRategy ",
06973 "EIGenvalue",
06974 "PAGe throw",
06975 "MINos errs",
06976 "EPSmachine",
06977 "OUTputfile",
06978 "BATch ",
06979 "INTeractiv",
06980 "VERsion ",
06981 "reserve ",
06982 "NODebug ",
06983 "DEBug ",
06984 "SHOw ",
06985 "SET "};
06986
06987 static Int_t nname = 25;
06988 static Int_t nntot = 30;
06989 static TString cprlev[5] = {
06990 "-1: NO OUTPUT EXCEPT FROM SHOW ",
06991 " 0: REDUCED OUTPUT ",
06992 " 1: NORMAL OUTPUT ",
06993 " 2: EXTRA OUTPUT FOR PROBLEM CASES",
06994 " 3: MAXIMUM OUTPUT "};
06995
06996 static TString cstrat[3] = {
06997 " 0: MINIMIZE THE NUMBER OF CALLS TO FUNCTION",
06998 " 1: TRY TO BALANCE SPEED AGAINST RELIABILITY",
06999 " 2: MAKE SURE MINIMUM TRUE, ERRORS CORRECT "};
07000
07001 static TString cdbopt[7] = {
07002 "REPORT ALL EXCEPTIONAL CONDITIONS ",
07003 "MNLINE: LINE SEARCH MINIMIZATION ",
07004 "MNDERI: FIRST DERIVATIVE CALCULATIONS ",
07005 "MNHESS: SECOND DERIVATIVE CALCULATIONS ",
07006 "MNMIGR: COVARIANCE MATRIX UPDATES ",
07007 "MNHES1: FIRST DERIVATIVE UNCERTAINTIES ",
07008 "MNCONT: MNCONTOUR PLOT (MNCROS SEARCH) "};
07009
07010
07011 Int_t f_inqu();
07012
07013
07014 Double_t val;
07015 Int_t iset, iprm, i, jseed, kname, iseed, iunit, id, ii, kk;
07016 Int_t ikseed, idbopt, igrain=0, iswsav, isw2;
07017 TString cfname, cmode, ckind, cwarn, copt, ctemp, ctemp2;
07018 Bool_t lname=kFALSE;
07019
07020 for (i = 1; i <= nntot; ++i) {
07021 ctemp = cname[i-1];
07022 ckind = ctemp(0,3);
07023 ctemp2 = fCword(4,6);
07024 if (strstr(ctemp2.Data(),ckind.Data())) goto L5;
07025 }
07026 i = 0;
07027 L5:
07028 kname = i;
07029
07030
07031 ctemp2 = fCword(0,3);
07032 if ( ctemp2.Contains("HEL")) goto L2000;
07033 if ( ctemp2.Contains("SHO")) goto L1000;
07034 if (!ctemp2.Contains("SET")) goto L1900;
07035
07036 ckind = "SET ";
07037
07038 if (kname <= 0) goto L1900;
07039
07040 switch ((int)kname) {
07041 case 1: goto L3000;
07042 case 2: goto L20;
07043 case 3: goto L30;
07044 case 4: goto L40;
07045 case 5: goto L3000;
07046 case 6: goto L60;
07047 case 7: goto L70;
07048 case 8: goto L80;
07049 case 9: goto L90;
07050 case 10: goto L100;
07051 case 11: goto L110;
07052 case 12: goto L120;
07053 case 13: goto L130;
07054 case 14: goto L140;
07055 case 15: goto L150;
07056 case 16: goto L160;
07057 case 17: goto L170;
07058 case 18: goto L3000;
07059 case 19: goto L190;
07060 case 20: goto L3000;
07061 case 21: goto L210;
07062 case 22: goto L220;
07063 case 23: goto L230;
07064 case 24: goto L240;
07065 case 25: goto L3000;
07066 case 26: goto L1900;
07067 case 27: goto L270;
07068 case 28: goto L280;
07069 case 29: goto L290;
07070 case 30: goto L300;
07071 }
07072
07073
07074 L20:
07075 iprm = Int_t(fWord7[0]);
07076 if (iprm > fNu) goto L25;
07077 if (iprm <= 0) goto L25;
07078 if (fNvarl[iprm-1] < 0) goto L25;
07079 fU[iprm-1] = fWord7[1];
07080 mnexin(fX);
07081 isw2 = fISW[1];
07082 mnrset(1);
07083
07084 fISW[1] = TMath::Min(isw2,1);
07085 fCfrom = "SET PARM";
07086 fNfcnfr = fNfcn;
07087 fCstatu = "NEW VALUES";
07088 return;
07089 L25:
07090 Printf(" UNDEFINED PARAMETER NUMBER. IGNORED.");
07091 return;
07092
07093 L30:
07094 mnlims();
07095 return;
07096
07097 L40:
07098
07099 goto L3000;
07100
07101 L60:
07102 fISW[4] = Int_t(fWord7[0]);
07103 return;
07104
07105 L70:
07106 fISW[2] = 0;
07107 return;
07108
07109 L80:
07110 mngrad();
07111 return;
07112
07113 L90:
07114 if (fWord7[0] == fUp) return;
07115 if (fWord7[0] <= 0) {
07116 if (fUp == fUpdflt) return;
07117 fUp = fUpdflt;
07118 } else {
07119 fUp = fWord7[0];
07120 }
07121 for (i = 1; i <= fNpar; ++i) {
07122 fErn[i-1] = 0;
07123 fErp[i-1] = 0;
07124 }
07125 mnwerr();
07126 return;
07127
07128
07129
07130 L100:
07131 goto L3000;
07132
07133 L110:
07134 fNpagwd = Int_t(fWord7[0]);
07135 fNpagwd = TMath::Max(fNpagwd,50);
07136 return;
07137
07138 L120:
07139 fNpagln = Int_t(fWord7[0]);
07140 return;
07141
07142
07143 L130:
07144 fLwarn = kFALSE;
07145 return;
07146
07147 L140:
07148 fLwarn = kTRUE;
07149 mnwarn("W", "SHO", "SHO");
07150 return;
07151
07152 L150:
07153 jseed = Int_t(fWord7[0]);
07154 val = 3;
07155 mnrn15(val, jseed);
07156 if (fISW[4] > 0) {
07157 Printf(" MINUIT RANDOM NUMBER SEED SET TO %d",jseed);
07158 }
07159 return;
07160
07161 L160:
07162
07163 goto L3000;
07164
07165 L170:
07166 fIstrat = Int_t(fWord7[0]);
07167 fIstrat = TMath::Max(fIstrat,0);
07168 fIstrat = TMath::Min(fIstrat,2);
07169 if (fISW[4] > 0) goto L1172;
07170 return;
07171
07172 L190:
07173 fNewpag = Int_t(fWord7[0]);
07174 goto L1190;
07175
07176 L210:
07177 if (fWord7[0] > 0 && fWord7[0] < .1) {
07178 fEpsmac = fWord7[0];
07179 }
07180 fEpsma2 = TMath::Sqrt(fEpsmac);
07181 goto L1210;
07182
07183 L220:
07184 iunit = Int_t(fWord7[0]);
07185 fIsyswr = iunit;
07186 fIstkwr[0] = iunit;
07187 if (fISW[4] >= 0) goto L1220;
07188 return;
07189
07190 L230:
07191 fISW[5] = 0;
07192 if (fISW[4] >= 0) goto L1100;
07193 return;
07194
07195 L240:
07196 fISW[5] = 1;
07197 if (fISW[4] >= 0) goto L1100;
07198 return;
07199
07200 L270:
07201 iset = 0;
07202 goto L281;
07203
07204 L280:
07205 iset = 1;
07206 L281:
07207 idbopt = Int_t(fWord7[0]);
07208 if (idbopt > 6) goto L288;
07209 if (idbopt >= 0) {
07210 fIdbg[idbopt] = iset;
07211 if (iset == 1) fIdbg[0] = 1;
07212 } else {
07213
07214 for (id = 0; id <= 6; ++id) { fIdbg[id] = iset; }
07215 }
07216 fLrepor = fIdbg[0] >= 1;
07217 mnwarn("D", "SHO", "SHO");
07218 return;
07219 L288:
07220 Printf(" UNKNOWN DEBUG OPTION %d REQUESTED. IGNORED",idbopt);
07221 return;
07222
07223 L290:
07224
07225 L300:
07226 goto L3000;
07227
07228 L1000:
07229
07230 ckind = "SHOW";
07231 if (kname <= 0) goto L1900;
07232
07233 switch ((int)kname) {
07234 case 1: goto L1010;
07235 case 2: goto L1020;
07236 case 3: goto L1030;
07237 case 4: goto L1040;
07238 case 5: goto L1050;
07239 case 6: goto L1060;
07240 case 7: goto L1070;
07241 case 8: goto L1070;
07242 case 9: goto L1090;
07243 case 10: goto L1100;
07244 case 11: goto L1110;
07245 case 12: goto L1120;
07246 case 13: goto L1130;
07247 case 14: goto L1130;
07248 case 15: goto L1150;
07249 case 16: goto L1160;
07250 case 17: goto L1170;
07251 case 18: goto L1180;
07252 case 19: goto L1190;
07253 case 20: goto L1200;
07254 case 21: goto L1210;
07255 case 22: goto L1220;
07256 case 23: goto L1100;
07257 case 24: goto L1100;
07258 case 25: goto L1250;
07259 case 26: goto L1900;
07260 case 27: goto L1270;
07261 case 28: goto L1270;
07262 case 29: goto L1290;
07263 case 30: goto L1300;
07264 }
07265
07266
07267 L1010:
07268 if (fAmin == fUndefi) mnamin();
07269 mnprin(0, fAmin);
07270 return;
07271
07272 L1020:
07273 if (fAmin == fUndefi) mnamin();
07274 mnprin(5, fAmin);
07275 return;
07276
07277 L1030:
07278 if (fAmin == fUndefi) mnamin();
07279 mnprin(1, fAmin);
07280 return;
07281
07282 L1040:
07283 mnmatu(1);
07284 return;
07285
07286 L1050:
07287 mnmatu(0);
07288 return;
07289
07290 L1060:
07291 if (fISW[4] < -1) fISW[4] = -1;
07292 if (fISW[4] > 3) fISW[4] = 3;
07293 Printf(" ALLOWED PRINT LEVELS ARE:");
07294 Printf(" %s",cprlev[0].Data());
07295 Printf(" %s",cprlev[1].Data());
07296 Printf(" %s",cprlev[2].Data());
07297 Printf(" %s",cprlev[3].Data());
07298 Printf(" %s",cprlev[4].Data());
07299 Printf(" CURRENT PRINTOUT LEVEL IS %s",cprlev[fISW[4]+1].Data());
07300 return;
07301
07302 L1070:
07303 if (fISW[2] <= 0) {
07304 Printf(" NOGRAD IS SET. DERIVATIVES NOT COMPUTED IN FCN.");
07305 } else {
07306 Printf(" GRAD IS SET. USER COMPUTES DERIVATIVES IN FCN.");
07307 }
07308 return;
07309
07310 L1090:
07311 Printf(" ERRORS CORRESPOND TO FUNCTION CHANGE OF %g",fUp);
07312 return;
07313
07314
07315 L1100:
07316
07317
07318
07319
07320
07321
07322
07323
07324
07325
07326
07327
07328
07329
07330
07331
07332
07333
07334
07335 cmode = "BATCH MODE ";
07336 if (fISW[5] == 1) cmode = "INTERACTIVE MODE";
07337 if (! lname) cfname = "unknown";
07338 Printf(" INPUT NOW BEING READ IN %s FROM UNIT NO. %d FILENAME: %s"
07339 ,(const char*)cmode,fIsysrd,(const char*)cfname);
07340 return;
07341
07342 L1110:
07343 Printf(" PAGE WIDTH IS SET TO %d COLUMNS",fNpagwd);
07344 return;
07345
07346 L1120:
07347 Printf(" PAGE LENGTH IS SET TO %d LINES",fNpagln);
07348 return;
07349
07350 L1130:
07351 cwarn = "SUPPRESSED";
07352 if (fLwarn) cwarn = "REPORTED ";
07353 Printf("%s",(const char*)cwarn);
07354 if (! fLwarn) mnwarn("W", "SHO", "SHO");
07355 return;
07356
07357 L1150:
07358 val = 0;
07359 mnrn15(val, igrain);
07360 ikseed = igrain;
07361 Printf(" MINUIT RNDM SEED IS CURRENTLY=%d",ikseed);
07362 val = 3;
07363 iseed = ikseed;
07364 mnrn15(val, iseed);
07365 return;
07366
07367 L1160:
07368 Printf(" TITLE OF CURRENT TASK IS:%s",(const char*)fCtitl);
07369 return;
07370
07371 L1170:
07372 Printf(" ALLOWED STRATEGIES ARE:");
07373 Printf(" %s",cstrat[0].Data());
07374 Printf(" %s",cstrat[1].Data());
07375 Printf(" %s",cstrat[2].Data());
07376 L1172:
07377 Printf(" NOW USING STRATEGY %s",(const char*)cstrat[fIstrat]);
07378 return;
07379
07380 L1180:
07381 iswsav = fISW[4];
07382 fISW[4] = 3;
07383 if (fISW[1] < 1) {
07384 Printf("%s",(const char*)fCovmes[0]);
07385 } else {
07386 mnpsdf();
07387 }
07388 fISW[4] = iswsav;
07389 return;
07390
07391 L1190:
07392 Printf(" PAGE THROW CARRIAGE CONTROL = %d",fNewpag);
07393 if (fNewpag == 0) {
07394 Printf(" NO PAGE THROWS IN MINUIT OUTPUT");
07395 }
07396 return;
07397
07398 L1200:
07399 for (ii = 1; ii <= fNpar; ++ii) {
07400 if (fErp[ii-1] > 0 || fErn[ii-1] < 0) goto L1204;
07401 }
07402 Printf(" THERE ARE NO MINOS ERRORS CURRENTLY VALID.");
07403 return;
07404 L1204:
07405 mnprin(4, fAmin);
07406 return;
07407
07408 L1210:
07409 Printf(" FLOATING-POINT NUMBERS ASSUMED ACCURATE TO %g",fEpsmac);
07410 return;
07411
07412 L1220:
07413 Printf(" MINUIT PRIMARY OUTPUT TO UNIT %d",fIsyswr);
07414 return;
07415
07416 L1250:
07417 Printf(" THIS IS MINUIT VERSION:%s",(const char*)fCvrsn);
07418 return;
07419
07420 L1270:
07421 for (id = 0; id <= 6; ++id) {
07422 copt = "OFF";
07423 if (fIdbg[id] >= 1) copt = "ON ";
07424 Printf(" DEBUG OPTION %3d IS %3s :%s"
07425 ,id,(const char*)copt,(const char*)cdbopt[id]);
07426 }
07427 if (! fLrepor) mnwarn("D", "SHO", "SHO");
07428 return;
07429
07430 L1290:
07431 ckind = "SHOW";
07432 goto L2100;
07433
07434 L1300:
07435 ckind = "SET ";
07436 goto L2100;
07437
07438
07439 L1900:
07440 Printf(" THE COMMAND:%10s IS UNKNOWN.",(const char*)fCword);
07441 goto L2100;
07442
07443
07444 L2000:
07445 ckind = "SET ";
07446 ctemp2 = fCword(3,7);
07447 if (strcmp(ctemp2.Data(), "SHO")) ckind = "SHOW";
07448 L2100:
07449 Printf(" THE FORMAT OF THE %4s COMMAND IS:",(const char*)ckind);
07450 Printf(" %s xxx [numerical arguments if any]",(const char*)ckind);
07451 Printf(" WHERE xxx MAY BE ONE OF THE FOLLOWING:");
07452 for (kk = 1; kk <= nname; ++kk) {
07453 Printf(" %s",cname[kk-1]);
07454 }
07455 return;
07456
07457
07458 L3000:
07459 Printf(" ABOVE COMMAND IS ILLEGAL. IGNORED");
07460
07461 }
07462
07463
07464 void TMinuit::mnsimp()
07465 {
07466
07467
07468
07469
07470
07471
07472
07473
07474 static Double_t alpha = 1;
07475 static Double_t beta = .5;
07476 static Double_t gamma = 2;
07477 static Double_t rhomin = 4;
07478 static Double_t rhomax = 8;
07479
07480
07481 Double_t dmin_, dxdi, yrho, f, ynpp1, aming, ypbar;
07482 Double_t bestx, ystar, y1, y2, ystst, pb, wg;
07483 Double_t absmin, rho, sig2, rho1, rho2;
07484 Int_t npfn, i, j, k, jhold, ncycl, nparx;
07485 Int_t nparp1, kg, jh, nf, jl, ns;
07486
07487 if (fNpar <= 0) return;
07488 if (fAmin == fUndefi) mnamin();
07489 fCfrom = "SIMPLEX ";
07490 fNfcnfr = fNfcn;
07491 fCstatu = "UNCHANGED ";
07492 npfn = fNfcn;
07493 nparp1 = fNpar + 1;
07494 nparx = fNpar;
07495 rho1 = alpha + 1;
07496 rho2 = rho1 + alpha*gamma;
07497 wg = 1 / Double_t(fNpar);
07498 if (fISW[4] >= 0) {
07499 Printf(" START SIMPLEX MINIMIZATION. CONVERGENCE WHEN EDM .LT. %g",fEpsi);
07500 }
07501 for (i = 1; i <= fNpar; ++i) {
07502 fDirin[i-1] = fWerr[i-1];
07503 mndxdi(fX[i-1], i-1, dxdi);
07504 if (dxdi != 0) fDirin[i-1] = fWerr[i-1] / dxdi;
07505 dmin_ = fEpsma2*TMath::Abs(fX[i-1]);
07506 if (fDirin[i-1] < dmin_) fDirin[i-1] = dmin_;
07507 }
07508
07509 L1:
07510 ynpp1 = fAmin;
07511 jl = nparp1;
07512 fSIMPy[nparp1-1] = fAmin;
07513 absmin = fAmin;
07514 for (i = 1; i <= fNpar; ++i) {
07515 aming = fAmin;
07516 fPbar[i-1] = fX[i-1];
07517 bestx = fX[i-1];
07518 kg = 0;
07519 ns = 0;
07520 nf = 0;
07521 L4:
07522 fX[i-1] = bestx + fDirin[i-1];
07523 mninex(fX);
07524 Eval(nparx, fGin, f, fU, 4); ++fNfcn;
07525 if (f <= aming) goto L6;
07526
07527 if (kg == 1) goto L8;
07528 kg = -1;
07529 ++nf;
07530 fDirin[i-1] *= -.4;
07531 if (nf < 3) goto L4;
07532 ns = 6;
07533
07534 L6:
07535 bestx = fX[i-1];
07536 fDirin[i-1] *= 3;
07537 aming = f;
07538 fCstatu = "PROGRESS ";
07539 kg = 1;
07540 ++ns;
07541 if (ns < 6) goto L4;
07542
07543 L8:
07544 fSIMPy[i-1] = aming;
07545 if (aming < absmin) jl = i;
07546 if (aming < absmin) absmin = aming;
07547 fX[i-1] = bestx;
07548 for (k = 1; k <= fNpar; ++k) { fP[k + i*fMaxpar - fMaxpar-1] = fX[k-1]; }
07549 }
07550 jh = nparp1;
07551 fAmin = fSIMPy[jl-1];
07552 mnrazz(ynpp1, fPbar, fSIMPy, jh, jl);
07553 for (i = 1; i <= fNpar; ++i) { fX[i-1] = fP[i + jl*fMaxpar - fMaxpar-1]; }
07554 mninex(fX);
07555 fCstatu = "PROGRESS ";
07556 if (fISW[4] >= 1) mnprin(5, fAmin);
07557 fEDM = fBigedm;
07558 sig2 = fEDM;
07559 ncycl = 0;
07560
07561 L50:
07562 if (sig2 < fEpsi && fEDM < fEpsi) goto L76;
07563 sig2 = fEDM;
07564 if (fNfcn - npfn > fNfcnmx) goto L78;
07565
07566 for (i = 1; i <= fNpar; ++i) {
07567 pb = 0;
07568 for (j = 1; j <= nparp1; ++j) { pb += wg*fP[i + j*fMaxpar - fMaxpar-1]; }
07569 fPbar[i-1] = pb - wg*fP[i + jh*fMaxpar - fMaxpar-1];
07570 fPstar[i-1] = (alpha + 1)*fPbar[i-1] - alpha*fP[i + jh*fMaxpar - fMaxpar-1];
07571 }
07572 mninex(fPstar);
07573 Eval(nparx, fGin, ystar, fU, 4); ++fNfcn;
07574 if (ystar >= fAmin) goto L70;
07575
07576 for (i = 1; i <= fNpar; ++i) {
07577 fPstst[i-1] = gamma*fPstar[i-1] + (1 - gamma)*fPbar[i-1];
07578 }
07579 mninex(fPstst);
07580 Eval(nparx, fGin, ystst, fU, 4); ++fNfcn;
07581
07582 y1 = (ystar - fSIMPy[jh-1])*rho2;
07583 y2 = (ystst - fSIMPy[jh-1])*rho1;
07584 rho = (rho2*y1 - rho1*y2)*.5 / (y1 - y2);
07585 if (rho < rhomin) goto L66;
07586 if (rho > rhomax) rho = rhomax;
07587 for (i = 1; i <= fNpar; ++i) {
07588 fPrho[i-1] = rho*fPbar[i-1] + (1 - rho)*fP[i + jh*fMaxpar - fMaxpar-1];
07589 }
07590 mninex(fPrho);
07591 Eval(nparx, fGin, yrho, fU, 4); ++fNfcn;
07592 if (yrho < fSIMPy[jl-1] && yrho < ystst) goto L65;
07593 if (ystst < fSIMPy[jl-1]) goto L67;
07594 if (yrho > fSIMPy[jl-1]) goto L66;
07595
07596 L65:
07597 mnrazz(yrho, fPrho, fSIMPy, jh, jl);
07598 goto L68;
07599 L66:
07600 if (ystst < fSIMPy[jl-1]) goto L67;
07601 mnrazz(ystar, fPstar, fSIMPy, jh, jl);
07602 goto L68;
07603 L67:
07604 mnrazz(ystst, fPstst, fSIMPy, jh, jl);
07605 L68:
07606 ++ncycl;
07607 if (fISW[4] < 2) goto L50;
07608 if (fISW[4] >= 3 || ncycl % 10 == 0) {
07609 mnprin(5, fAmin);
07610 }
07611 goto L50;
07612
07613 L70:
07614 if (ystar >= fSIMPy[jh-1]) goto L73;
07615 jhold = jh;
07616 mnrazz(ystar, fPstar, fSIMPy, jh, jl);
07617 if (jhold != jh) goto L50;
07618
07619 L73:
07620 for (i = 1; i <= fNpar; ++i) {
07621 fPstst[i-1] = beta*fP[i + jh*fMaxpar - fMaxpar-1] + (1 - beta)*fPbar[i-1];
07622 }
07623 mninex(fPstst);
07624 Eval(nparx, fGin, ystst, fU, 4); ++fNfcn;
07625 if (ystst > fSIMPy[jh-1]) goto L1;
07626
07627 if (ystst < fAmin) goto L67;
07628 mnrazz(ystst, fPstst, fSIMPy, jh, jl);
07629 goto L50;
07630
07631 L76:
07632 if (fISW[4] >= 0) {
07633 Printf(" SIMPLEX MINIMIZATION HAS CONVERGED.");
07634 }
07635 fISW[3] = 1;
07636 goto L80;
07637 L78:
07638 if (fISW[4] >= 0) {
07639 Printf(" SIMPLEX TERMINATES WITHOUT CONVERGENCE.");
07640 }
07641 fCstatu = "CALL LIMIT";
07642 fISW[3] = -1;
07643 fISW[0] = 1;
07644 L80:
07645 for (i = 1; i <= fNpar; ++i) {
07646 pb = 0;
07647 for (j = 1; j <= nparp1; ++j) { pb += wg*fP[i + j*fMaxpar - fMaxpar-1]; }
07648 fPbar[i-1] = pb - wg*fP[i + jh*fMaxpar - fMaxpar-1];
07649 }
07650 mninex(fPbar);
07651 Eval(nparx, fGin, ypbar, fU, 4); ++fNfcn;
07652 if (ypbar < fAmin) mnrazz(ypbar, fPbar, fSIMPy, jh, jl);
07653 mninex(fX);
07654 if (fNfcnmx + npfn - fNfcn < fNpar*3) goto L90;
07655 if (fEDM > fEpsi*2) goto L1;
07656 L90:
07657 if (fISW[4] >= 0) mnprin(5, fAmin);
07658 }
07659
07660
07661 void TMinuit::mnstat(Double_t &fmin, Double_t &fedm, Double_t &errdef, Int_t &npari, Int_t &nparx, Int_t &istat)
07662 {
07663
07664
07665
07666
07667
07668
07669
07670
07671
07672
07673
07674
07675
07676
07677
07678
07679 fmin = fAmin;
07680 fedm = fEDM;
07681 errdef = fUp;
07682 npari = fNpar;
07683 nparx = fNu;
07684 istat = fISW[1];
07685 if (fEDM == fBigedm) fedm = fUp;
07686 if (fAmin == fUndefi) {
07687 fmin = 0;
07688 fedm = fUp;
07689 istat = 0;
07690 }
07691 }
07692
07693
07694 void TMinuit::mntiny(Double_t epsp1, Double_t &epsbak)
07695 {
07696
07697
07698
07699
07700
07701
07702
07703 epsbak = epsp1 - 1;
07704 }
07705
07706
07707 Bool_t TMinuit::mnunpt(TString &cfname)
07708 {
07709
07710
07711
07712
07713 Int_t i, l, ic;
07714 Bool_t ret_val;
07715 static TString cpt = " ABCDEFGHIJKLMNOPQRSTUVWXYZabcdefghijklmnopqrstuvwxyz1234567890./;:[]$%*_!@#&+()";
07716
07717 ret_val = kFALSE;
07718 l = strlen((const char*)cfname);
07719 for (i = 1; i <= l; ++i) {
07720 for (ic = 1; ic <= 80; ++ic) {
07721 if (cfname[i-1] == cpt[ic-1]) goto L100;
07722 }
07723 return kTRUE;
07724 L100:
07725 ;
07726 }
07727 return ret_val;
07728 }
07729
07730
07731 void TMinuit::mnvert(Double_t *a, Int_t l, Int_t, Int_t n, Int_t &ifail)
07732 {
07733
07734
07735
07736
07737
07738
07739
07740
07741 Int_t a_offset;
07742
07743
07744 Double_t si;
07745 Int_t i, j, k, kp1, km1;
07746
07747
07748 a_offset = l + 1;
07749 a -= a_offset;
07750
07751
07752 ifail = 0;
07753 if (n < 1) goto L100;
07754 if (n > fMaxint) goto L100;
07755
07756 for (i = 1; i <= n; ++i) {
07757 si = a[i + i*l];
07758 if (si <= 0) goto L100;
07759 fVERTs[i-1] = 1 / TMath::Sqrt(si);
07760 }
07761 for (i = 1; i <= n; ++i) {
07762 for (j = 1; j <= n; ++j) {
07763 a[i + j*l] = a[i + j*l]*fVERTs[i-1]*fVERTs[j-1];
07764 }
07765 }
07766
07767 for (i = 1; i <= n; ++i) {
07768 k = i;
07769
07770 if (a[k + k*l] != 0) fVERTq[k-1] = 1 / a[k + k*l];
07771 else goto L100;
07772 fVERTpp[k-1] = 1;
07773 a[k + k*l] = 0;
07774 kp1 = k + 1;
07775 km1 = k - 1;
07776 if (km1 < 0) goto L100;
07777 else if (km1 == 0) goto L50;
07778 else goto L40;
07779 L40:
07780 for (j = 1; j <= km1; ++j) {
07781 fVERTpp[j-1] = a[j + k*l];
07782 fVERTq[j-1] = a[j + k*l]*fVERTq[k-1];
07783 a[j + k*l] = 0;
07784 }
07785 L50:
07786 if (k - n < 0) goto L51;
07787 else if (k - n == 0) goto L60;
07788 else goto L100;
07789 L51:
07790 for (j = kp1; j <= n; ++j) {
07791 fVERTpp[j-1] = a[k + j*l];
07792 fVERTq[j-1] = -a[k + j*l]*fVERTq[k-1];
07793 a[k + j*l] = 0;
07794 }
07795
07796 L60:
07797 for (j = 1; j <= n; ++j) {
07798 for (k = j; k <= n; ++k) { a[j + k*l] += fVERTpp[j-1]*fVERTq[k-1]; }
07799 }
07800 }
07801
07802 for (j = 1; j <= n; ++j) {
07803 for (k = 1; k <= j; ++k) {
07804 a[k + j*l] = a[k + j*l]*fVERTs[k-1]*fVERTs[j-1];
07805 a[j + k*l] = a[k + j*l];
07806 }
07807 }
07808 return;
07809
07810 L100:
07811 ifail = 1;
07812 }
07813
07814
07815 void TMinuit::mnwarn(const char *copt1, const char *corg1, const char *cmes1)
07816 {
07817
07818
07819
07820
07821
07822
07823
07824
07825
07826
07827
07828
07829 TString copt = copt1;
07830 TString corg = corg1;
07831 TString cmes = cmes1;
07832
07833 const Int_t kMAXMES = 10;
07834 Int_t ityp, i, ic, nm;
07835 TString englsh, ctyp;
07836
07837 if (corg(0,3) != "SHO" || cmes(0,3) != "SHO") {
07838
07839
07840 if (copt == "W") {
07841 ityp = 1;
07842 if (fLwarn) {
07843 Printf(" MINUIT WARNING IN %s",(const char*)corg);
07844 Printf(" ============== %s",(const char*)cmes);
07845 return;
07846 }
07847 } else {
07848 ityp = 2;
07849 if (fLrepor) {
07850 Printf(" MINUIT DEBUG FOR %s",(const char*)corg);
07851 Printf(" =============== %s ",(const char*)cmes);
07852 return;
07853 }
07854 }
07855
07856 if (fNwrmes[ityp-1] == 0) fIcirc[ityp-1] = 0;
07857 ++fNwrmes[ityp-1];
07858 ++fIcirc[ityp-1];
07859 if (fIcirc[ityp-1] > 10) fIcirc[ityp-1] = 1;
07860 ic = fIcirc[ityp-1];
07861 fOrigin[ic] = corg;
07862 fWarmes[ic] = cmes;
07863 fNfcwar[ic] = fNfcn;
07864 return;
07865 }
07866
07867
07868 if (copt == "W") {
07869 ityp = 1;
07870 ctyp = "WARNING";
07871 } else {
07872 ityp = 2;
07873 ctyp = "*DEBUG*";
07874 }
07875 if (fNwrmes[ityp-1] > 0) {
07876 englsh = " WAS SUPPRESSED. ";
07877 if (fNwrmes[ityp-1] > 1) englsh = "S WERE SUPPRESSED.";
07878 Printf(" %5d MINUIT %s MESSAGE%s",fNwrmes[ityp-1]
07879 ,(const char*)ctyp,(const char*)englsh);
07880 nm = fNwrmes[ityp-1];
07881 ic = 0;
07882 if (nm > kMAXMES) {
07883 Printf(" ONLY THE MOST RECENT 10 WILL BE LISTED BELOW.");
07884 nm = kMAXMES;
07885 ic = fIcirc[ityp-1];
07886 }
07887 Printf(" CALLS ORIGIN MESSAGE");
07888 for (i = 1; i <= nm; ++i) {
07889 ++ic;
07890 if (ic > kMAXMES) ic = 1;
07891 Printf(" %6d %s %s", fNfcwar[ic],fOrigin[ic].Data(),fWarmes[ic].Data());
07892 }
07893 fNwrmes[ityp-1] = 0;
07894 Printf(" ");
07895 }
07896 }
07897
07898
07899 void TMinuit::mnwerr()
07900 {
07901
07902
07903
07904
07905
07906
07907 Double_t denom, ba, al, dx, du1, du2;
07908 Int_t ndex, ierr, i, j, k, l, ndiag, k1, iin;
07909
07910
07911 if (fISW[1] >= 1) {
07912 for (l = 1; l <= fNpar; ++l) {
07913 ndex = l*(l + 1) / 2;
07914 dx = TMath::Sqrt(TMath::Abs(fVhmat[ndex-1]*fUp));
07915 i = fNexofi[l-1];
07916 if (fNvarl[i-1] > 1) {
07917 al = fAlim[i-1];
07918 ba = fBlim[i-1] - al;
07919 du1 = al + 0.5*(TMath::Sin(fX[l-1] + dx) + 1)*ba - fU[i-1];
07920 du2 = al + 0.5*(TMath::Sin(fX[l-1] - dx) + 1)*ba - fU[i-1];
07921 if (dx > 1) du1 = ba;
07922 dx = 0.5*(TMath::Abs(du1) + TMath::Abs(du2));
07923 }
07924 fWerr[l-1] = dx;
07925 }
07926 }
07927
07928 if (fISW[1] >= 1) {
07929 for (i = 1; i <= fNpar; ++i) {
07930 fGlobcc[i-1] = 0;
07931 k1 = i*(i-1) / 2;
07932 for (j = 1; j <= i; ++j) {
07933 k = k1 + j;
07934 fP[i + j*fMaxpar - fMaxpar-1] = fVhmat[k-1];
07935 fP[j + i*fMaxpar - fMaxpar-1] = fP[i + j*fMaxpar - fMaxpar-1];
07936 }
07937 }
07938 mnvert(fP, fMaxint, fMaxint, fNpar, ierr);
07939 if (ierr == 0) {
07940 for (iin = 1; iin <= fNpar; ++iin) {
07941 ndiag = iin*(iin + 1) / 2;
07942 denom = fP[iin + iin*fMaxpar - fMaxpar-1]*fVhmat[ndiag-1];
07943 if (denom <= 1 && denom >= 0) fGlobcc[iin-1] = 0;
07944 else fGlobcc[iin-1] = TMath::Sqrt(1 - 1 / denom);
07945 }
07946 }
07947 }
07948 }