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