35 const char charal[29] =
" .ABCDEFGHIJKLMNOPQRSTUVWXYZ";
62 m_userParameterIdToInternalId(2*maxpar+1),
63 m_userParameterValue(2*maxpar+1,0),
65 m_userParameterFlag(2*maxpar+1,-1),
249 for (
int i = 0; i <
fMaxpar; i++) {
294 const string& sname = name;
295 mnparm( parNo, sname, initVal, initErr, lowerLimit, upperLimit, err);
414 if (
fFCN) (*fFCN)(npar,grad,fval,par,flag);
423 Int_urt kint = m_userParameterIdToInternalId[externalParameterId];
441 mnexcm(
"FIX", &tmp, 1, err );
455 mnpout( parNo, name, currentValue, currentError, bnd1, bnd2, err );
495 mnexcm(
"MIGRAD", plist, npar, err );
510 mnexcm(
"HESSE", plist, npar, err );
525 mnexcm(
"MINOS", plist, npar, err );
539 mnexcm(
"RELEASE", &tmp, 1, err );
549 mnexcm(
"SET ERRDEF", &up, 1, err );
569 mnexcm(
"SET PRINT", &tmp, 1, err );
590 *m_logStream <<
" FIRST CALL TO USER FUNCTION AT NEW START POINT, WITH IFLAG=4." << endl;;
613 Int_urt kwid, lwid, na=0, log_;
617 if (al == ah) ah = al + 1;
620 if (naa == -1)
goto L150;
629 if (awid <= 1) --log_;
632 if (sigfig > 2)
goto L40;
636 if (sigfig > 2.5)
goto L50;
640 if (sigfig > 5)
goto L60;
651 if (bwid <= 0)
goto L10;
662 if (naa > 5)
goto L240;
663 if (naa == -1)
return;
665 if (naa > 1 || nb == 1)
return;
670 if (nb << 1 != naa)
return;
686 Int_urt ndex, i, j, m, n, nparx;
693 for (i = 1; i <=
fNpar; ++i) {
695 for (j = 1; j <=
fNpar; ++j) {
698 ndex = m*(m-1) / 2 + n;
703 for (i = 1; i <=
fNpar; ++i) {denom +=
fGrd[i-1]*(
fXt[i-1] - pvec[i-1]); }
709 ycalf = (f -
fApsi) / denom;
728 for (i = 1; i <=
fMaxext; ++i) {
730 m_userParameterValue[i] = 0;
732 m_userParameterName[i] =
fCundef;
734 m_userParameterFlag[i] = -1;
736 m_userParameterIdToInternalId[i] = 0;
755 static string clabel =
"0123456789ABCDEFGHIJ";
760 Double_urt ylabel, fmn, fmx, xlo, ylo, xup, yup;
761 Double_urt devs, xsav, ysav, bwidx, bwidy, unext, ff, xb4;
762 Int_urt i, ngrid, ixmid, nparx, ix, nx, ny, ki1, ki2, ixzero, iy, ics;
763 string chmid, chln, chzero;
768 if (ke1 <= 0 || ke2 <= 0)
goto L1350;
769 if (ke1 >
fNu || ke2 >
fNu)
goto L1350;
772 ki1 = m_userParameterIdToInternalId[ke1];
773 ki2 = m_userParameterIdToInternalId[ke2];
774 if (ki1 <= 0 || ki2 <= 0)
goto L1350;
775 if (ki1 == ki2)
goto L1350;
784 xsav = m_userParameterValue[ke1];
785 ysav = m_userParameterValue[ke2];
787 if (devs <= 0) devs = 2;
792 xlo = m_userParameterValue[ke1] - devs*
fWerr[ki1-1];
793 xup = m_userParameterValue[ke1] + devs*fWerr[ki1-1];
794 ylo = m_userParameterValue[ke2] - devs*fWerr[ki2-1];
795 yup = m_userParameterValue[ke2] + devs*fWerr[ki2-1];
807 if (nx < 11) nx = 11;
808 if (ny < 11) ny = 11;
809 if (nx >= 115) nx = 114;
813 if (m_userParameterFlag[ke1] > 1) {
818 if (m_userParameterFlag[ke2] > 1) {
827 for (i = 1; i <= 20; ++i) { contur[i-1] =
fAmin +
fUp*(i-1)*(i-1); }
828 contur[0] +=
fUp*.01;
831 m_userParameterValue[ke2] = yup;
835 chmid.resize(nx+1,
' ');
836 chzero.resize(nx+1,
' ');
837 chln.resize(nx+1,
' ');
838 for (ix = 1; ix <= nx + 1; ++ix) {
840 m_userParameterValue[ke1] = xlo +
Double_urt(ix-1)*bwidx;
842 Eval(nparx,
fGin, ff, m_userParameterValue, 4);
845 if (xb4 < 0 && m_userParameterValue[ke1] > 0) ixzero = ix - 1;
847 xb4 = m_userParameterValue[ke1];
853 *m_logStream <<
" Y-AXIS: PARAMETER " << setw(3) << ke2 <<
": " << m_userParameterName[ke2] <<
'\n';
855 chzero[ixzero-1] =
'+';
857 *m_logStream <<
" X=0\n";
861 for (iy = 1; iy <= ny; ++iy) {
863 unext = m_userParameterValue[ke2] - bwidy;
867 chln.resize(nx+1,
' ');
869 if (ixzero != 0) chln[ixzero-1] =
':';
872 if (m_userParameterValue[ke2] > ysav && unext < ysav) chln = chmid;
873 if (m_userParameterValue[ke2] > 0 && unext < 0) chln = chzero;
876 m_userParameterValue[ke2] = unext;
877 ylabel = m_userParameterValue[ke2] + bwidy*.5;
879 for (ix = 1; ix <= nx + 1; ++ix) {
880 fcna[ix-1] = fcnb[ix-1];
882 m_userParameterValue[ke1] = xlo +
Double_urt(ix-1)*bwidx;
884 Eval(nparx,
fGin, ff, m_userParameterValue, 4);
888 for (ix = 1; ix <= nx; ++ix) {
895 for (ics = 1; ics <= 20; ++ics) {
896 if (contur[ics-1] > fmn)
goto L240;
900 if (contur[ics-1] < fmx) chln[ix-1] = clabel[ics-1];
903 *m_logStream << setw(12) << setprecision(4) << ylabel << chln <<
'\n';
908 chln.resize(nx+1,
' ');
909 chln.replace(0,1,
"I");
910 chln.replace(ixmid-1,1,
"I");
911 chln.replace(nx-1,1,
"I");
913 *m_logStream <<
" " << chln <<
'\n';
919 *m_logStream <<
" " << setw(12) << setprecision(4) << xlo << chln << setw(12) << xup <<
'\n';
921 *m_logStream <<
" " << chln << setw(12) << setprecision(4) << xsav <<
'\n';
924 *m_logStream <<
" " << setw(12) << setprecision(4) << xlo << chln
925 << setw(12) << setprecision(4) << xsav << chln
926 << setw(12) << setprecision(4) << xup <<
'\n';
932 *m_logStream <<
" X-AXIS: PARAMETER" << setw(3) << ke1 <<
": "<< m_userParameterName[ke1]
933 <<
" ONE COLUMN=" << setw(12) << setprecision(4) << bwidx <<
'\n';
935 *m_logStream <<
" FUNCTION VALUES: F(I)=" << setw(12) << setprecision(4) <<
fAmin 936 <<
" +" << setw(12) << setprecision(4) <<
fUp <<
" *I**2\n";
940 m_userParameterValue[ke1] = xsav;
941 m_userParameterValue[ke2] = ysav;
943 *m_logStream << endl;
946 *m_logStream <<
" INVALID PARAMETER NUMBER(S) REQUESTED. IGNORED." << endl;
977 Int_urt ierr, ipos, i, llist, lenbuf, lnc;
979 string comand, crdbuf, ctemp;
982 int (*pf)(int)=toupper;
983 transform(crdbuf.begin(), crdbuf.end(), crdbuf.begin(), pf);
984 lenbuf = crdbuf.length();
990 if (crdbuf[i-1] ==
'\'')
break;
991 if (crdbuf[i-1] ==
' ') {
1000 *m_logStream <<
" BLANK COMMAND IGNORED." << endl;
1007 if (crdbuf.substr(ipos-1,3) ==
"PAR") {
1013 if (crdbuf.substr(ipos-1,7) ==
"SET INP") {
1019 if (crdbuf.substr(ipos-1,7) ==
"SET TIT") {
1025 if (crdbuf.substr(ipos-1,7) ==
"SET COV") {
1031 ctemp = crdbuf.substr(ipos-1,lenbuf-ipos+1);
1035 *m_logStream <<
" COMMAND CANNOT BE INTERPRETED" << endl;
1066 Double_urt dist, xdir, ydir, aopt, u1min, u2min;
1068 Double_urt a1, a2, val2mi, val2pl, dc, sclfac, bigdis, sigsav;
1069 Int_urt nall, iold, line, mpar, ierr, inew, move, next, i, j, nfcol, iercr;
1070 Int_urt idist=0, npcol, kints, i2, i1, lr, nfcnco=0, ki1, ki2, ki3, ke3;
1071 Int_urt nowpts, istrav, nfmxin, isw2, isw4;
1079 ldebug =
fIdbg[6] >= 1;
1082 if (ke1 <= 0 || ke2 <= 0)
goto L1350;
1083 if (ke1 >
fNu || ke2 >
fNu)
goto L1350;
1086 ki1 = m_userParameterIdToInternalId[ke1];
1087 ki2 = m_userParameterIdToInternalId[ke2];
1088 if (ki1 <= 0 || ki2 <= 0)
goto L1350;
1089 if (ki1 == ki2)
goto L1350;
1090 if (nptu < 4)
goto L1400;
1098 u1min = m_userParameterValue[ke1];
1099 u2min = m_userParameterValue[ke2];
1105 *m_logStream <<
" START MNCONTOUR CALCULATION OF" << setw(4) << nptu <<
" POINTS ON CONTOUR." << endl;
1108 ki3 = 6 - ki1 - ki2;
1112 *m_logStream <<
" EACH POINT IS A MINIMUM WITH RESPECT TO PARAMETER " << setw(3) << ke3
1113 <<
" " << m_userParameterName[ke3] << endl;
1116 *m_logStream <<
" EACH POINT IS A MINIMUM WITH RESPECT TO THE OTHER" << setw(3) << (
fNpar - 2)
1117 <<
" VARIABLE PARAMETERS." << endl;
1124 mnmnot(ke1, ke2, val2pl, val2mi);
1126 xptu[0] =
fAlim[ke1-1];
1127 mnwarn(
"W",
"MNContour ",
"Contour squeezed by parameter limits.");
1129 if (
fErn[ki1-1] >= 0)
goto L1500;
1130 xptu[0] = u1min +
fErn[ki1-1];
1135 xptu[2] =
fBlim[ke1-1];
1136 mnwarn(
"W",
"MNContour ",
"Contour squeezed by parameter limits.");
1138 if (
fErp[ki1-1] <= 0)
goto L1500;
1139 xptu[2] = u1min +
fErp[ki1-1];
1142 scalx = 1 / (xptu[2] - xptu[0]);
1144 mnmnot(ke2, ke1, val2pl, val2mi);
1146 yptu[1] =
fAlim[ke2-1];
1147 mnwarn(
"W",
"MNContour ",
"Contour squeezed by parameter limits.");
1149 if (
fErn[ki2-1] >= 0)
goto L1500;
1150 yptu[1] = u2min +
fErn[ki2-1];
1154 yptu[3] =
fBlim[ke2-1];
1155 mnwarn(
"W",
"MNContour ",
"Contour squeezed by parameter limits.");
1157 if (
fErp[ki2-1] <= 0)
goto L1500;
1158 yptu[3] = u2min +
fErp[ki2-1];
1161 scaly = 1 / (yptu[3] - yptu[1]);
1166 *m_logStream <<
" Plot of four points found by MINOS" << endl;
1172 for (i = 2; i <= nall; ++i) {
1173 fXpt[i-1] = xptu[i-2];
1174 fYpt[i-1] = yptu[i-2];
1177 *m_logStream <<
fChpt <<
" ABCD" << endl;
1191 for (i = 1; i <= mpar; ++i) {
fXt[i-1] =
fX[i-1]; }
1192 i__1 = mpar*(mpar + 1) / 2;
1194 for (i = 1; i <= mpar; ++i) {
1200 kints = m_userParameterIdToInternalId[ke1];
1203 kints = m_userParameterIdToInternalId[ke2];
1206 for (inew = next; inew <= nptu; ++inew) {
1209 for (iold = 1; iold <= inew - 1; ++iold) {
1211 if (i2 == inew) i2 = 1;
1212 d__1 = scalx*(xptu[iold-1] - xptu[i2-1]);
1213 d__2 = scaly*(yptu[iold-1] - yptu[i2-1]);
1214 dist = d__1*d__1 + d__2*d__2;
1215 if (dist > bigdis) {
1222 if (i2 == inew) i2 = 1;
1227 fXmidcr = a1*xptu[i1-1] + a2*xptu[i2-1];
1228 fYmidcr = a1*yptu[i1-1] + a2*yptu[i2-1];
1229 xdir = yptu[i2-1] - yptu[i1-1];
1230 ydir = xptu[i1-1] - xptu[i2-1];
1244 *m_logStream <<
" MNCONT CANNOT FIND NEXT POINT ON CONTOUR. ONLY" << setw(3) << nowpts
1245 <<
" POINTS FOUND." << endl;
1249 mnwarn(
"W",
"MNContour ",
"Cannot find midpoint, try closer.");
1255 for (move = nowpts; move >= i1 + 1; --move) {
1256 xptu[move] = xptu[move-1];
1257 yptu[move] = yptu[move-1];
1267 if (nowpts < nptu)
fCstatu =
"INCOMPLETE";
1275 for (i = 2; i <= nall; ++i) {
1276 fXpt[i-1] = xptu[i-2];
1277 fYpt[i-1] = yptu[i-2];
1283 *m_logStream <<
" Y-AXIS: PARAMETER " << setw(3) << ke2 <<
" " << m_userParameterName[ke2] <<
'\n';
1289 *m_logStream <<
" X-AXIS: PARAMETER " << setw(3) << ke1 <<
" " << m_userParameterName[ke1] <<
'\n';
1293 npcol = (nowpts + 1) / 2;
1296 *m_logStream << setw(5) << nowpts <<
" POINTS ON CONTOUR. FMIN=" << setw(13) << setprecision(5) << abest
1297 <<
" ERRDEF=" << setw(11) << setprecision(3) <<
fUp <<
'\n';
1306 *m_logStream <<
" " << m_userParameterName[ke1] << m_userParameterName[ke2]
1307 << m_userParameterName[ke1] << m_userParameterName[ke1] <<
'\n';
1308 for (line = 1; line <= nfcol; ++line) {
1311 *m_logStream << setw(6) << line << setw(13) << setprecision(5) << xptu[line-1] << setw(13) << setprecision(5) << yptu[line-1]
1312 << setw(15) << lr << setw(13) << setprecision(5) << xptu[lr-1] << setw(13) << setprecision(5) << yptu[lr-1] <<
'\n';
1314 if (nfcol < npcol) {
1316 *m_logStream << setw(6) << npcol << setw(13) << setprecision(5) << xptu[npcol-1] << setw(13) << setprecision(5) << yptu[npcol-1] <<
'\n';
1323 i__1 = mpar*(mpar + 1) / 2;
1325 for (i = 1; i <= mpar; ++i) {
1341 m_userParameterValue[ke1] = u1min;
1342 m_userParameterValue[ke2] = u2min;
1347 *m_logStream <<
" INVALID PARAMETER NUMBERS.\n";
1351 *m_logStream <<
" LESS THAN FOUR POINTS REQUESTED.\n";
1358 *m_logStream <<
" MNCONT UNABLE TO FIND FOUR POINTS.\n";
1361 m_userParameterValue[ke1] = u1min;
1362 m_userParameterValue[ke2] = u2min;
1368 *m_logStream << endl;
1393 const char *cnumer =
"123456789-.0+";
1396 Int_urt ifld, iend, lend, left, nreq, ipos, kcmnd, nextb, ic, ibegin, ltoadd;
1397 Int_urt ielmnt, lelmnt[25], nelmnt;
1403 char *crdbuf =
const_cast<char*
>( cardbuf.c_str() );
1404 lend = cardbuf.length();
1410 for (ipos = nextb; ipos <= lend; ++ipos) {
1412 if (crdbuf[ipos-1] ==
' ')
continue;
1413 if (crdbuf[ipos-1] ==
',')
goto L250;
1419 for (ipos = ibegin + 1; ipos <= lend; ++ipos) {
1420 if (crdbuf[ipos-1] ==
' ')
goto L250;
1421 if (crdbuf[ipos-1] ==
',')
goto L250;
1427 if (iend >= ibegin) celmnt[ielmnt-1] = &crdbuf[ibegin-1];
1428 else celmnt[ielmnt-1] = cnull;
1429 lelmnt[ielmnt-1] = iend - ibegin + 1;
1430 if (lelmnt[ielmnt-1] > 19) {
1432 *m_logStream <<
" MINUIT WARNING: INPUT DATA WORD TOO LONG." << endl;
1433 ctemp = cardbuf.substr(ibegin-1,iend-ibegin+1);
1435 *m_logStream <<
" ORIGINAL:" << ctemp <<
'\n';
1437 *m_logStream <<
" TRUNCATED TO:" << celmnt[ielmnt-1] << endl;
1438 lelmnt[ielmnt-1] = 19;
1440 if (ipos >= lend)
goto L300;
1441 if (ielmnt >= 25)
goto L300;
1443 for (ipos = iend + 1; ipos <= lend; ++ipos) {
1444 if (crdbuf[ipos-1] ==
' ')
continue;
1446 if (crdbuf[ipos-1] ==
',') nextb = ipos + 1;
1453 command[0] =
' '; command[1] = 0;
1457 if (ielmnt == 0)
goto L900;
1459 for (ielmnt = 1; ielmnt <= nelmnt; ++ielmnt) {
1460 if ( celmnt[ielmnt-1] == cnull)
goto L450;
1461 for (ic = 1; ic <= 13; ++ic) {
1462 if (*celmnt[ielmnt-1] == cnumer[ic-1])
goto L450;
1464 if (kcmnd >= maxcwd)
continue;
1465 left = maxcwd - kcmnd;
1466 ltoadd = lelmnt[ielmnt-1];
1467 if (ltoadd > left) ltoadd = left;
1468 strncpy(&command[kcmnd],celmnt[ielmnt-1],ltoadd);
1470 if (kcmnd == maxcwd)
continue;
1471 command[kcmnd] =
' ';
1481 for (ifld = ielmnt; ifld <= nelmnt; ++ifld) {
1484 nreq = nelmnt - ielmnt + 1;
1486 *m_logStream <<
" MINUIT WARNING IN MNCRCK: \n";
1488 *m_logStream <<
" COMMAND HAS INPUT" << setw(5) << nreq <<
" NUMERIC FIELDS, BUT MINUIT CAN ACCEPT ONLY" 1489 << setw(3) << mxp << endl;
1492 if (celmnt[ifld-1] == cnull) plist[llist-1] = 0;
1494 sscanf(celmnt[ifld-1],
"%lf",&plist[llist-1]);
1499 if (lnc <= 0) lnc = 1;
1518 Double_urt alsb[3], flsb[3], bmin, bmax, zmid, sdev, zdir, zlim;
1519 Double_urt coeff[3], aleft, aulim, fdist, adist, aminsv;
1520 Double_urt anext, fnext, slope, s1, s2, x1, x2, ecarmn, ecarmx;
1521 Double_urt determ, rt, smalla, aright, aim, tla, tlf, dfda,ecart;
1522 Int_urt iout=0, i, ileft, ierev, maxlk, ibest, ik, it;
1523 Int_urt noless, iworst=0, iright, itoohi, kex, ipt;
1528 ldebug =
fIdbg[6] >= 1;
1547 for (ik = 1; ik <= 2; ++ik) {
1553 if (
fKe2cr == 0)
continue;
1559 if (m_userParameterFlag[kex] <= 1)
continue;
1560 if (zdir == 0)
continue;
1561 zlim =
fAlim[kex-1];
1562 if (zdir > 0) zlim =
fBlim[kex-1];
1571 mneval(anext, fnext, ierev);
1574 std::printf(
" MNCROS: calls=%8d AIM=%10.5f F,A=%10.5f%10.5f\n",
fNfcn,aim,fnext,aopt);
1576 if (ierev > 0)
goto L900;
1577 if (
fLimset && fnext <= aim)
goto L930;
1579 fXpt[ipt-1] = anext;
1580 fYpt[ipt-1] = fnext;
1588 if (aopt < -.5)aopt = -.5;
1589 if (aopt > 1) aopt = 1;
1595 mneval(aopt, fnext, ierev);
1598 std::printf(
" MNCROS: calls=%8d AIM=%10.5f F,A=%10.5f%10.5f\n",
fNfcn,aim,fnext,aopt);
1600 if (ierev > 0)
goto L900;
1601 if (
fLimset && fnext <= aim)
goto L930;
1604 fXpt[ipt-1] = alsb[1];
1605 fYpt[ipt-1] = fnext;
1608 dfda = (flsb[1] - flsb[0]) / (alsb[1] - alsb[0]);
1610 if (dfda > 0)
goto L460;
1612 mnwarn(
"D",
"MNCROS ",
"Looking for slope of the right sign");
1614 for (it = 1; it <= maxlk; ++it) {
1623 mneval(aopt, fnext, ierev);
1626 std::printf(
" MNCROS: calls=%8d AIM=%10.5f F,A=%10.5f%10.5f\n",
fNfcn,aim,fnext,aopt);
1628 if (ierev > 0)
goto L900;
1629 if (
fLimset && fnext <= aim)
goto L930;
1632 fXpt[ipt-1] = alsb[1];
1633 fYpt[ipt-1] = fnext;
1636 dfda = (flsb[1] - flsb[0]) / (alsb[1] - alsb[0]);
1637 if (dfda > 0)
goto L450;
1639 mnwarn(
"W",
"MNCROS ",
"Cannot find slope of the right sign");
1644 aopt = alsb[1] + (aim - flsb[1]) / dfda;
1649 if (adist < tla && fdist < tlf)
goto L800;
1650 if (ipt >= 15)
goto L950;
1652 if (aopt < bmin) aopt = bmin;
1654 if (aopt > bmax) aopt = bmax;
1661 mneval(aopt, fnext, ierev);
1664 std::printf(
" MNCROS: calls=%8d AIM=%10.5f F,A=%10.5f%10.5f\n",
fNfcn,aim,fnext,aopt);
1666 if (ierev > 0)
goto L900;
1667 if (
fLimset && fnext <= aim)
goto L930;
1670 fXpt[ipt-1] = alsb[2];
1671 fYpt[ipt-1] = fnext;
1679 for (i = 1; i <= 3; ++i) {
1681 if (ecart > ecarmx) { ecarmx = ecart; iworst = i; }
1682 if (ecart < ecarmn) { ecarmn = ecart; ibest = i; }
1683 if (flsb[i-1] < aim) ++noless;
1686 if (noless == 1 || noless == 2)
goto L500;
1688 if (noless == 0 && ibest != 3)
goto L950;
1691 if (noless == 3 && ibest != 3) {
1697 alsb[iworst-1] = alsb[2];
1698 flsb[iworst-1] = flsb[2];
1699 dfda = (flsb[1] - flsb[0]) / (alsb[1] - alsb[0]);
1703 mnpfit(alsb, flsb, 3, coeff, sdev);
1704 if (coeff[2] <= 0) {
1705 mnwarn(
"D",
"MNCROS ",
"Curvature is negative near contour line.");
1707 determ = coeff[1]*coeff[1] - coeff[2]*4*(coeff[0] - aim);
1709 mnwarn(
"D",
"MNCROS ",
"Problem 2, impossible determinant");
1714 x1 = (-coeff[1] + rt) / (coeff[2]*2);
1715 x2 = (-coeff[1] - rt) / (coeff[2]*2);
1716 s1 = coeff[1] + x1*2*coeff[2];
1717 s2 = coeff[1] + x2*2*coeff[2];
1720 *m_logStream <<
" MNCONTour problem 1" << endl;
1734 if (ipt >= 15)
goto L950;
1743 for (i = 1; i <= 3; ++i) {
1745 if (ecart < ecarmn) { ecarmn = ecart; ibest = i; }
1746 if (ecart > ecarmx) { ecarmx = ecart; }
1747 if (flsb[i-1] > aim) {
1748 if (iright == 0) iright = i;
1749 else if (flsb[i-1] > flsb[iright-1]) iout = i;
1750 else { iout = iright; iright = i; }
1752 else if (ileft == 0) ileft = i;
1753 else if (flsb[i-1] < flsb[ileft-1]) iout = i;
1754 else { iout = ileft; ileft = i; }
1757 if (ecarmx >
URMath::Abs(flsb[iout-1] - aim)*10) {
1758 aopt = aopt*.5 + (alsb[iright-1] + alsb[ileft-1])*.25;
1762 if (slope*smalla > tlf) smalla = tlf / slope;
1763 aleft = alsb[ileft-1] + smalla;
1764 aright = alsb[iright-1] - smalla;
1766 if (aopt < aleft) aopt = aleft;
1767 if (aopt > aright) aopt = aright;
1768 if (aleft > aright) aopt = (aleft + aright)*.5;
1777 mneval(aopt, fnext, ierev);
1780 std::printf(
" MNCROS: calls=%8d AIM=%10.5f F,A=%10.5f%10.5f\n",
fNfcn,aim,fnext,aopt);
1782 if (ierev > 0)
goto L900;
1783 if (
fLimset && fnext <= aim)
goto L930;
1786 fYpt[ipt-1] = fnext;
1789 alsb[iout-1] = aopt;
1790 flsb[iout-1] = fnext;
1802 if (ierev == 1)
goto L940;
1819 for (i = 1; i <= ipt; ++i) {
1820 if (
fYpt[i-1] > aim + fUp) {
1828 if (
fXdircr < 0) chsign =
"NEGA";
1831 *m_logStream << chsign <<
"TIVE MINOS ERROR, PARAMETER " << setw(3) <<
fKe1cr << endl;
1835 *m_logStream <<
"POINTS LABELLED '+' WERE TOO HIGH TO PLOT." << endl;
1839 *m_logStream <<
"RIGHTMOST POINT IS UP AGAINST LIMIT." << endl;
1861 *m_logStream <<
" FUNCTION MUST BE MINIMIZED BEFORE CALLING " <<
fCfrom << endl;
1868 mnwarn(
"W",
fCfrom.c_str(),
"NO ERROR MATRIX. WILL IMPROVISE.\n");
1869 for (i = 1; i <=
fNpar; ++i) {
1871 for (j = 1; j <= i-1; ++j) {
1876 if (
fG2[i-1] <= 0) {
1880 if (m_userParameterFlag[iext] > 1) {
1885 fG2[i-1] =
fUp / (wint*wint);
1908 Double_urt tlrstp, tlrgrd, epspri, optstp, stpmax, stpmin, fs2, grbfor=0, d1d2, xtf;
1909 Int_urt icyc, ncyc, iint, iext, i, nparx;
1913 ldebug =
fIdbg[2] >= 1;
1915 if (
fISW[2] == 1)
goto L100;
1925 ostringstream warning;
1926 warning <<
"function value differs from AMIN by " << df;
1927 mnwarn(
"D",
"MNDERI", warning.str().c_str());
1932 *m_logStream <<
" FIRST DERIVATIVE DEBUG PRINTOUT. MNDERI\n";
1933 *m_logStream <<
" PAR DERIV STEP MINSTEP OPTSTEP D1-D2 2ND DRV" << endl;
1950 for (i = 1; i <=
fNpar; ++i) {
1957 for (icyc = 1; icyc <= ncyc; ++icyc) {
1963 if (
fGstep[i-1] < 0 && step > .5) step = .5;
1966 if (step > stpmax) step = stpmax;
1969 if (step < stpmin) step = stpmin;
1971 if (URMath::Abs((step - stepb4) / step) < tlrstp)
goto L50;
1977 fX[i-1] = xtf + step;
1982 fX[i-1] = xtf - step;
1987 fGrd[i-1] = (fs1 - fs2) / (step*2);
1988 fG2[i-1] = (fs1 + fs2 -
fAmin*2) / (step*step);
1991 d1d2 = (fs1 + fs2 -
fAmin*2) / step;
1992 std::printf(
"%4d%11.3g%11.3g%10.2g%10.2g%10.2g%10.2g\n",i,
fGrd[i-1],step,stpmin,optstp,d1d2,
fG2[i-1]);
1995 if (URMath::Abs(grbfor -
fGrd[i-1]) / (URMath::Abs(
fGrd[i-1]) + dfmin/step) < tlrgrd)
1999 if (ncyc == 1)
goto L50;
2000 { ostringstream warning2;
2001 warning2 <<
"First derivative not converged. " <<
fGrd[i-1] << grbfor;
2002 mnwarn(
"D",
"MNDERI", warning2.str().c_str()); }
2010 for (iint = 1; iint <=
fNpar; ++iint) {
2013 if (m_userParameterFlag[iext] <= 1) {
2035 if (m_userParameterFlag[i] > 1) {
2051 Int_urt i, j, k, l, m=0, i0, i1, j1, m1, n1;
2055 a_offset = ndima + 1;
2063 for (i1 = 2; i1 <= n; ++i1) {
2065 f = a[i + (i-1)*ndima];
2068 if (l < 1)
goto L25;
2070 for (k = 1; k <= l; ++k) {
2071 d__1 = a[i + k*ndima];
2077 if (gl > 1e-35)
goto L30;
2085 if (f >= 0) gl = -gl;
2088 a[i + (i-1)*ndima] = f - gl;
2090 for (j = 1; j <= l; ++j) {
2091 a[j + i*ndima] = a[i + j*ndima] / h;
2093 for (k = 1; k <= j; ++k) { gl += a[j + k*ndima]*a[i + k*ndima]; }
2094 if (j >= l)
goto L47;
2096 for (k = j1; k <= l; ++k) { gl += a[k + j*ndima]*a[i + k*ndima]; }
2098 work[n + j] = gl / h;
2099 f += gl*a[j + i*ndima];
2102 for (j = 1; j <= l; ++j) {
2104 gl = work[n + j] - hh*f;
2106 for (k = 1; k <= j; ++k) {
2107 a[j + k*ndima] = a[j + k*ndima] - f*work[n + k] - gl*a[i + k*ndima];
2116 for (i = 1; i <= n; ++i) {
2118 if (work[i] == 0 || l == 0)
goto L100;
2120 for (j = 1; j <= l; ++j) {
2122 for (k = 1; k <= l; ++k) { gl += a[i + k*ndima]*a[k + j*ndima]; }
2123 for (k = 1; k <= l; ++k) { a[k + j*ndima] -= gl*a[k + i*ndima]; }
2126 work[i] = a[i + i*ndima];
2128 if (l == 0)
continue;
2130 for (j = 1; j <= l; ++j) {
2137 for (i = 2; i <= n; ++i) {
2139 work[i0] = work[i0 + 1];
2144 for (l = 1; l <= n; ++l) {
2148 for (m1 = l; m1 <= n; ++m1) {
2154 if (m == l)
goto L205;
2157 if (j == mits)
return;
2159 pt = (work[l + 1] - work[l]) / (work[n + l]*2);
2162 if (pt < 0) pr = pt - r;
2164 h = work[l] - work[n + l] / pr;
2165 for (i = l; i <= n; ++i) { work[i] -= h; }
2172 for (i1 = l; i1 <= m1; ++i1) {
2179 c = pt / work[n + i];
2181 work[n + j] = s*work[n + i]*r;
2186 c = work[n + i] / pt;
2188 work[n + j] = s*pt*r;
2192 pt = c*work[i] - s*gl;
2193 work[j] = h + s*(c*gl + s*work[i]);
2194 for (k = 1; k <= n; ++k) {
2196 a[k + j*ndima] = s*a[k + i*ndima] + c*h;
2197 a[k + i*ndima] = c*a[k + i*ndima] - s*h;
2208 for (i = 1; i <= n1; ++i) {
2212 for (j = i1; j <= n; ++j) {
2213 if (work[j] >= pt)
continue;
2218 if (k == i)
continue;
2222 for (j = 1; j <= n; ++j) {
2223 pt = a[j + i*ndima];
2224 a[j + i*ndima] = a[j + k*ndima];
2225 a[j + k*ndima] = pt;
2241 Int_urt emat_dim1, emat_offset;
2245 Int_urt i, j, k, npard, k2, kk, iz, nperln, kga, kgb;
2250 emat_offset = emat_dim1 + 1;
2251 emat -= emat_offset;
2254 if (
fISW[1] < 1)
return;
2257 *m_logStream <<
" EXTERNAL ERROR MATRIX. NDIM=" << setw(4) << ndim <<
" NPAR=" << setw(3)
2258 <<
fNpar <<
" ERR DEF=" <<
fUp << endl;
2265 std::printf(
" USER-DIMENSIONED ARRAY EMAT NOT BIG ENOUGH. REDUCED MATRIX CALCULATED.\n");
2272 if (
fISW[4] >= 1 && npard > nperln) {
2278 for (i = 1; i <= npard; ++i) {
2281 for (j = 1; j <= i; ++j) {
2284 emat[i + j*emat_dim1] = dxdi*
fVhmat[kgb-1]*dxdj*
fUp;
2285 emat[j + i*emat_dim1] = emat[i + j*emat_dim1];
2290 for (i = 1; i <= npard; ++i) {
2292 if (npard >= nperln) iz = i;
2293 for (k = 1; nperln < 0 ? k >= iz : k <= iz; k += nperln) {
2294 k2 = k + nperln - 1;
2295 if (k2 > iz) k2 = iz;
2297 for (kk = k; kk <= k2; ++kk) {
2298 ostringstream localTemp;
2299 localTemp <<
" " << setw(10) << setprecision(3) << emat[i + kk*emat_dim1];
2300 ctemp += localTemp.str();
2302 std::printf(
"%s\n",ctemp.c_str());
2328 if (iex >
fNu || iex <= 0)
goto L900;
2331 iin = m_userParameterIdToInternalId[iex];
2333 if (iin <= 0)
goto L900;
2336 eplus =
fErp[iin-1];
2337 if (eplus ==
fUndefi) eplus = 0;
2338 eminus =
fErn[iin-1];
2339 if (eminus ==
fUndefi) eminus = 0;
2341 ndiag = iin*(iin + 1) / 2;
2345 if (
fISW[1] < 2)
return;
2387 if (
fISW[0] >= 1) ierev = 1;
2388 if (
fISW[3] < 1) ierev = 2;
2416 string comand = command;
2417 static string clower =
"abcdefghijklmnopqrstuvwxyz";
2418 static string cupper =
"ABCDEFGHIJKLMNOPQRSTUVWXYZ";
2419 const string cname[40] = {
2464 Double_urt step, xptu[101], yptu[101], f, rno;
2465 Int_urt icol, kcol, ierr, iint, iext, lnow, nptu, i, iflag, ierrf;
2466 Int_urt ilist, nparx, izero, nf, lk, it, iw, inonde, nsuper;
2467 Int_urt it2, ke1, ke2, nowprt, kll, krl;
2468 string chwhy, c26, cvblnk, cneway, comd;
2476 lk = comand.length();
2477 if (lk > 20) lk = 20;
2479 int (*pf)(int)=toupper;
2483 for (iw = 1; iw <=
fMaxpar; ++iw) {
2485 if (iw <= llist)
fWord7[iw-1] = plist[iw-1];
2489 if (
fCword.substr(0,7) !=
"SET PRI" ||
fWord7[0] >= 0) {
2492 if (lnow > 4) lnow = 4;
2493 std::printf(
" **********\n");
2494 ostringstream octemp;
2496 ctemp = octemp.str();
2498 for (i = 1; i <= lnow; ++i) {
2499 ostringstream localTemp;
2500 localTemp <<
" " << plist[i-1];
2501 ctemp += localTemp.str();
2503 std::printf(
"%s\n",ctemp.c_str());
2511 std::printf(
" ***********\n");
2512 for (i = lnow + 1; i <= kll; ++i) {
2513 std::printf(
"%12.4g\n",plist[i-1]);
2516 std::printf(
" **********\n");
2518 std::printf(
" ERROR: ABOVE CALL TO MNEXCM TRIED TO PASS MORE THAN %d PARAMETERS.\n",
fMaxpar);
2535 ctemp =
fCword.substr(0,3);
2536 for (i = 1; i <= nntot; ++i) {
2537 if (strncmp(ctemp.c_str(),cname[i-1].c_str(),3) == 0)
goto L90;
2539 std::printf(
"UNKNOWN COMMAND IGNORED:%s\n", comand.c_str());
2544 if (
fCword.substr(0,4) ==
"MINO") i = 5;
2545 if (i != 6 && i != 7 && i != 8 && i != 23) {
2560 case 10:
goto L1000;
2561 case 11:
goto L1100;
2562 case 12:
goto L1200;
2563 case 13:
goto L1300;
2564 case 14:
goto L1400;
2565 case 15:
goto L1500;
2566 case 16:
goto L1600;
2567 case 17:
goto L1700;
2568 case 18:
goto L1800;
2569 case 19:
goto L1900;
2570 case 20:
goto L1900;
2571 case 21:
goto L1900;
2572 case 22:
goto L2200;
2573 case 23:
goto L2300;
2574 case 24:
goto L2400;
2575 case 25:
goto L1900;
2576 case 26:
goto L2600;
2577 case 27:
goto L3300;
2578 case 28:
goto L3300;
2579 case 29:
goto L3300;
2580 case 30:
goto L3300;
2581 case 31:
goto L3300;
2582 case 32:
goto L3300;
2583 case 33:
goto L3300;
2584 case 34:
goto L3400;
2585 case 35:
goto L3500;
2586 case 36:
goto L3600;
2587 case 37:
goto L3700;
2588 case 38:
goto L3800;
2589 case 39:
goto L3900;
2590 case 40:
goto L4000;
2599 if (
fISW[3] < 1) ierflg = 4;
2607 if (
fISW[3] >= 1)
return;
2609 if (
fISW[0] == 1)
return;
2610 if (
fCword.substr(0,3) ==
"MIG")
return;
2615 if (
fISW[0] == 1)
return;
2618 if (
fISW[3] >= 1) ierflg = 0;
2633 if (fNfcn < nsuper)
goto L510;
2634 std::printf(
" TOO MANY FUNCTION CALLS. MINOS GIVES UP\n");
2654 std::printf(
"%s: NO PARAMETERS REQUESTED \n",
fCword.c_str());
2657 for (ilist = 1; ilist <= llist; ++ilist) {
2658 iext =
Int_urt(plist[ilist-1]);
2659 chwhy =
" IS UNDEFINED.";
2660 if (iext <= 0)
goto L930;
2661 if (iext >
fNu)
goto L930;
2663 if (m_userParameterFlag[iext] < 0)
goto L930;
2664 chwhy =
" IS CONSTANT. ";
2666 if (m_userParameterFlag[iext] == 0)
goto L930;
2668 iint = m_userParameterIdToInternalId[iext];
2670 chwhy =
" ALREADY FIXED.";
2671 if (iint == 0)
goto L930;
2673 if (ierr == 0) lfixed =
kurTRUE;
2676 chwhy =
" ALREADY VARIABLE.";
2677 if (iint > 0)
goto L930;
2684 std::printf(
" PARAMETER%4d %s IGNORED.\n",iext,chwhy.c_str());
2686 if (lfreed || lfixed)
mnrset(0);
2699 if (it > 1 || it < 0)
goto L1005;
2710 std::printf(
" IGNORED. UNKNOWN ARGUMENT:%4d\n",it);
2720 if (iext <= 0)
goto L1210;
2723 if (iext <=
fNu) it2 = m_userParameterIdToInternalId[iext];
2724 if (it2 <= 0)
goto L1250;
2730 std::printf(
" PARAMETER%4d NOT VARIABLE.\n",iext);
2742 std::printf(
"%s: NO PARAMETERS REQUESTED \n",
fCword.c_str());
2750 if (ierrf > 0) ierflg = 3;
2782 }
else if (f <
fAmin) {
2786 if (
fISW[4] >= 0 && iflag <= 5 && nowprt == 1) {
2789 if (iflag == 3)
fFval3 = f;
2791 if (iflag > 5)
mnrset(1);
2802 std::printf(
" CALL TO USER FUNCTION WITH IFLAG = 3\n");
2808 if (
fCword.substr(0,3) ==
"END") ierflg = 10;
2809 if (
fCword.substr(0,3) ==
"RET") ierflg = 12;
2815 std::printf(
" MINUIT MEMORY CLEARED. NO PARAMETERS NOW DEFINED.\n");
2821 for (icol = 5; icol <= lk; ++icol) {
2822 if (
fCword[icol-1] ==
' ')
continue;
2827 if (kcol == 0) comd =
"* ";
2828 else comd =
fCword.substr(kcol-1,lk-kcol+1);
2836 if (ke1 == 0 &&
fNpar == 2) {
2841 if (nptu <= 0) nptu = 20;
2842 if (nptu > 101) nptu = 101;
2845 mncont(ke1, ke2, nptu, xptu, yptu, ierrf);
2846 if (ierrf < nptu) ierflg = 4;
2847 if (ierrf == -1) ierflg = 3;
2852 if (step <= 0) step = 2;
2855 for (i = 1; i <=
fNpar; ++i) {
2858 fX[i-1] += rno*step*
fWerr[i-1];
2866 std::printf(
" BLANK COMMAND IGNORED.\n");
2872 std::printf(
" THE *COVARIANCE* COMMAND IS OSBSOLETE. THE COVARIANCE MATRIX IS NOW SAVED IN A DIFFERENT FORMAT WITH THE *SAVE* COMMAND AND READ IN WITH:*SET COVARIANCE*\n");
2877 cneway =
"SET PRInt ";
2881 cneway =
"SET GRAd ";
2885 cneway =
"SHOW COVar";
2889 cneway =
"SET ERRdef";
2893 cneway =
"SET LIMits";
2900 std::printf(
" OBSOLETE COMMAND:%s PLEASE USE: %s\n",
fCword.c_str()
2903 if (fCword ==
"SAVE ")
goto L1500;
2921 for (iint = 1; iint <=
fNpar; ++iint) {
2924 mnpint(m_userParameterValue[iext], iext, pinti);
2925 pint[iint-1] = pinti;
2939 Int_urt kold, nold, ndex, knew, iext, i, j, m, n, lc, ik;
2944 if (iint >
fNpar || iint <= 0) {
2946 std::printf(
" MINUIT ERROR. ARGUMENT TO MNFIXP=%4d\n",iint);
2952 std::printf(
" MINUIT CANNOT FIX PARAMETER%4d MAXIMUM NUMBER THAT CAN BE FIXED IS %d\n",iext,
fMaxpar);
2958 m_userParameterIdToInternalId[iext] = 0;
2973 for (ik = iext + 1; ik <=
fNu; ++ik) {
2975 if (m_userParameterIdToInternalId[ik] > 0) {
2977 lc = m_userParameterIdToInternalId[ik] - 1;
2979 m_userParameterIdToInternalId[ik] = lc;
2990 if (
fISW[1] <= 0)
return;
2992 if (
fNpar <= 0)
return;
2993 for (i = 1; i <= nold; ++i) {
2996 ndex = m*(m-1) / 2 + n;
3002 for (i = 1; i <= nold; ++i) {
3003 for (j = 1; j <= i; ++j) {
3005 if (j == iint || i == iint)
continue;
3030 Double_urt grdv, xv, dirinv, g2v, gstepv, xtv;
3031 Int_urt i, ipsav, ka, lc, ik, iq, ir, is;
3034 std::printf(
" CALL TO MNFREE IGNORED. ARGUMENT GREATER THAN ONE\n");
3037 std::printf(
" CALL TO MNFREE IGNORED. THERE ARE NO FIXED PARAMETERS\n");
3039 if (k == 1 || k == 0)
goto L40;
3044 if (m_userParameterIdToInternalId[ka] == 0)
goto L15;
3045 std::printf(
" IGNORED. PARAMETER SPECIFIED IS ALREADY VARIABLE.\n");
3048 if (
fNpfix < 1)
goto L21;
3049 for (ik = 1; ik <=
fNpfix; ++ik) {
if (
fIpfix[ik-1] == ka)
goto L24; }
3051 std::printf(
" PARAMETER%4d NOT FIXED. CANNOT BE RELEASED.\n",ka);
3054 if (ik ==
fNpfix)
goto L40;
3064 for (i = ik + 1; i <=
fNpfix; ++i) {
3082 if (
fNpfix < 1)
goto L300;
3086 for (ik =
fNu; ik >= ir; --ik) {
3088 if (m_userParameterIdToInternalId[ik] > 0) {
3090 lc = m_userParameterIdToInternalId[ik] + 1;
3093 m_userParameterIdToInternalId[ik] = lc;
3095 fX[lc-1] =
fX[lc-2];
3105 if (is == 0) is =
fNpar;
3107 m_userParameterIdToInternalId[ir] = is;
3110 fX[is-1] =
fXs[iq-1];
3123 std::printf(
" PARAMETER%4d %s RESTORED TO VARIABLE.\n",ir,
3124 m_userParameterName[ir].c_str());
3126 if (k == 0)
goto L40;
3148 static string cwd =
" ";
3152 if (
fWord7[0] > 0)
goto L2000;
3167 std::printf(
" CHECK OF GRADIENT CALCULATION IN FCN\n");
3168 std::printf(
" PARAMETER G(IN FCN) G(MINUIT) DG(MINUIT) AGREEMENT\n");
3171 for (lc = 1; lc <=
fNpar; ++lc) {
3181 if (cwd !=
"GOOD")
fISW[2] = 0;
3185 std::printf(
" %5d %10s%12.4e%12.4e%12.4e %s\n",i
3186 ,(m_userParameterName[i]).c_str()
3190 std::printf(
" AGREEMENT=NONE MEANS FCN DID NOT CALCULATE THE DERIVATIVE\n");
3193 std::printf(
" MINUIT DOES NOT ACCEPT DERIVATIVE CALCULATIONS BY FCN\n");
3194 std::printf(
" TO FORCE ACCEPTANCE, ENTER *SET GRAD 1*\n");
3229 string comd = comd_in;
3230 int (*pf)(int)=toupper;
3231 transform(comd.begin(), comd.end(), comd.begin(), pf);
3233 if( comd.length() == 0 || comd.substr(0,1) ==
"*" || comd.substr(0,1) ==
"?" || comd==
"HELP" ) {
3234 std::printf(
" ==>List of MINUIT Interactive commands:\n");
3235 std::printf(
" CLEar Reset all parameter names and values undefined\n");
3236 std::printf(
" CONtour Make contour map of the user function\n");
3237 std::printf(
" EXIT Exit from Interactive Minuit\n");
3238 std::printf(
" FIX Cause parameter(s) to remain constant\n");
3239 std::printf(
" HESse Calculate the Hessian or error matrix.\n");
3240 std::printf(
" IMPROVE Search for a new minimum around current minimum\n");
3241 std::printf(
" MIGrad Minimize by the method of Migrad\n");
3242 std::printf(
" MINImize MIGRAD + SIMPLEX method if Migrad fails\n");
3243 std::printf(
" MINOs Exact (non-linear) parameter error analysis\n");
3244 std::printf(
" MNContour Calculate one MINOS function contour\n");
3245 std::printf(
" PARameter Define or redefine new parameters and values\n");
3246 std::printf(
" RELease Make previously FIXed parameters variable again\n");
3247 std::printf(
" REStore Release last parameter fixed\n");
3248 std::printf(
" SAVe Save current parameter values on a file\n");
3249 std::printf(
" SCAn Scan the user function by varying parameters\n");
3250 std::printf(
" SEEk Minimize by the method of Monte Carlo\n");
3251 std::printf(
" SET Set various MINUIT constants or conditions\n");
3252 std::printf(
" SHOw Show values of current constants or conditions\n");
3253 std::printf(
" SIMplex Minimize by the method of Simplex\n");
3262 if( !strncmp(comd.c_str(),
"CLE",3) ) {
3263 std::printf(
" ***>CLEAR");
3264 std::printf(
" Resets all parameter names and values to undefined.\n");
3265 std::printf(
" Must normally be followed by a PARameters command or \n");
3266 std::printf(
" equivalent, in order to define parameter values.\n");
3274 if( !strncmp(comd.c_str(),
"CON",3) ) {
3275 std::printf(
" ***>CONTOUR <par1> <par2> [devs] [ngrid]\n");
3276 std::printf(
" Instructs Minuit to trace contour lines of the user function\n");
3277 std::printf(
" with respect to the two parameters whose external numbers\n");
3278 std::printf(
" are <par1> and <par2>.\n");
3279 std::printf(
" Other variable parameters of the function, if any, will have\n");
3280 std::printf(
" their values fixed at the current values during the contour\n");
3281 std::printf(
" tracing. The optional parameter [devs] (default value 2.)\n");
3282 std::printf(
" gives the number of standard deviations in each parameter\n");
3283 std::printf(
" which should lie entirely within the plotting area.\n");
3284 std::printf(
" Optional parameter [ngrid] (default value 25 unless page\n");
3285 std::printf(
" size is too small) determines the resolution of the plot,\n");
3286 std::printf(
" i.e. the number of rows and columns of the grid at which the\n");
3287 std::printf(
" function will be evaluated. [See also MNContour.]\n");
3295 if( !strncmp(comd.c_str(),
"END",3) ) {
3296 std::printf(
" ***>END\n");
3297 std::printf(
" Signals the end of a data block (i.e., the end of a fit),\n");
3298 std::printf(
" and implies that execution should continue, because another\n");
3299 std::printf(
" Data Block follows. A Data Block is a set of Minuit data\n");
3300 std::printf(
" consisting of\n");
3301 std::printf(
" (1) A Title,\n");
3302 std::printf(
" (2) One or more Parameter Definitions,\n");
3303 std::printf(
" (3) A blank line, and\n");
3304 std::printf(
" (4) A set of Minuit Commands.\n");
3305 std::printf(
" The END command is used when more than one Data Block is to\n");
3306 std::printf(
" be used with the same FCN function. It first causes Minuit\n");
3307 std::printf(
" to issue a CALL FCN with IFLAG=3, in order to allow FCN to\n");
3308 std::printf(
" perform any calculations associated with the final fitted\n");
3309 std::printf(
" parameter values, unless a CALL FCN 3 command has already\n");
3310 std::printf(
" been executed at the current FCN value.\n");
3318 if( !strncmp(comd.c_str(),
"EXI",3) ) {
3319 std::printf(
" ***>EXIT\n");
3320 std::printf(
" Signals the end of execution.\n");
3321 std::printf(
" The EXIT command first causes Minuit to issue a CALL FCN\n");
3322 std::printf(
" with IFLAG=3, to allow FCN to perform any calculations\n");
3323 std::printf(
" associated with the final fitted parameter values, unless a\n");
3324 std::printf(
" CALL FCN 3 command has already been executed.\n");
3332 if( !strncmp(comd.c_str(),
"FIX",3) ) {
3333 std::printf(
" ***>FIX} <parno> [parno] ... [parno]\n");
3334 std::printf(
" Causes parameter(s) <parno> to be removed from the list of\n");
3335 std::printf(
" variable parameters, and their value(s) will remain constant\n");
3336 std::printf(
" during subsequent minimizations, etc., until another command\n");
3337 std::printf(
" changes their value(s) or status.\n");
3345 if( !strncmp(comd.c_str(),
"HES",3) ) {
3346 std::printf(
" ***>HESse [maxcalls]");
3347 std::printf(
" Calculate, by finite differences, the Hessian or error matrix.\n");
3348 std::printf(
" That is, it calculates the full matrix of second derivatives\n");
3349 std::printf(
" of the function with respect to the currently variable\n");
3350 std::printf(
" parameters, and inverts it, printing out the resulting error\n");
3351 std::printf(
" matrix. The optional argument [maxcalls] specifies the\n");
3352 std::printf(
" (approximate) maximum number of function calls after which\n");
3353 std::printf(
" the calculation will be stopped.\n");
3361 if( !strncmp(comd.c_str(),
"IMP",3) ) {
3362 std::printf(
" ***>IMPROVE [maxcalls]");
3363 std::printf(
" If a previous minimization has converged, and the current\n");
3364 std::printf(
" values of the parameters therefore correspond to a local\n");
3365 std::printf(
" minimum of the function, this command requests a search for\n");
3366 std::printf(
" additional distinct local minima.\n");
3367 std::printf(
" The optional argument [maxcalls] specifies the (approximate\n");
3368 std::printf(
" maximum number of function calls after which the calculation\n");
3369 std::printf(
" will be stopped.\n");
3377 if( !strncmp(comd.c_str(),
"MIG",3) ) {
3378 std::printf(
" ***>MIGrad [maxcalls] [tolerance]\n");
3379 std::printf(
" Causes minimization of the function by the method of Migrad,\n");
3380 std::printf(
" the most efficient and complete single method, recommended\n");
3381 std::printf(
" for general functions (see also MINImize).\n");
3382 std::printf(
" The minimization produces as a by-product the error matrix\n");
3383 std::printf(
" of the parameters, which is usually reliable unless warning\n");
3384 std::printf(
" messages are produced.\n");
3385 std::printf(
" The optional argument [maxcalls] specifies the (approximate)\n");
3386 std::printf(
" maximum number of function calls after which the calculation\n");
3387 std::printf(
" will be stopped even if it has not yet converged.\n");
3388 std::printf(
" The optional argument [tolerance] specifies required tolerance\n");
3389 std::printf(
" on the function value at the minimum.\n");
3390 std::printf(
" The default tolerance is 0.1, and the minimization will stop\n");
3391 std::printf(
" when the estimated vertical distance to the minimum (EDM) is\n");
3392 std::printf(
" less than 0.001*[tolerance]*UP (see [SET ERRordef]).\n");
3400 if( !strncmp(comd.c_str(),
"MINI",4) ) {
3401 std::printf(
" ***>MINImize [maxcalls] [tolerance]\n");
3402 std::printf(
" Causes minimization of the function by the method of Migrad,\n");
3403 std::printf(
" as does the MIGrad command, but switches to the SIMplex method\n");
3404 std::printf(
" if Migrad fails to converge. Arguments are as for MIGrad.\n");
3405 std::printf(
" Note that command requires four characters to be unambiguous.\n");
3413 if( !strncmp(comd.c_str(),
"MIN0",4) ) {
3414 std::printf(
" ***>MINOs [maxcalls] [parno] [parno] ...\n");
3415 std::printf(
" Causes a Minos error analysis to be performed on the parameters\n");
3416 std::printf(
" whose numbers [parno] are specified. If none are specified,\n");
3417 std::printf(
" Minos errors are calculated for all variable parameters.\n");
3418 std::printf(
" Minos errors may be expensive to calculate, but are very\n");
3419 std::printf(
" reliable since they take account of non-linearities in the\n");
3420 std::printf(
" problem as well as parameter correlations, and are in general\n");
3421 std::printf(
" asymmetric.\n");
3422 std::printf(
" The optional argument [maxcalls] specifies the (approximate)\n");
3423 std::printf(
" maximum number of function calls per parameter requested,\n");
3424 std::printf(
" after which the calculation will stop for that parameter.\n");
3432 if( !strncmp(comd.c_str(),
"MNC",3) ) {
3433 std::printf(
" ***>MNContour <par1> <par2> [npts]\n");
3434 std::printf(
" Calculates one function contour of FCN with respect to\n");
3435 std::printf(
" parameters par1 and par2, with FCN minimized always with\n");
3436 std::printf(
" respect to all other NPAR-2 variable parameters (if any).\n");
3437 std::printf(
" Minuit will try to find npts points on the contour (default 20)\n");
3438 std::printf(
" If only two parameters are variable at the time, it is not\n");
3439 std::printf(
" necessary to specify their numbers. To calculate more than\n");
3440 std::printf(
" one contour, it is necessary to SET ERRordef to the appropriate\n");
3441 std::printf(
" value and issue the MNContour command for each contour.\n");
3449 if( !strncmp(comd.c_str(),
"PAR",3) ) {
3450 std::printf(
" ***>PARameters\n");
3451 std::printf(
" followed by one or more parameter definitions.\n");
3452 std::printf(
" Parameter definitions are of the form:\n");
3453 std::printf(
" <number> ''name'' <value> <step> [lolim] [uplim] \n");
3454 std::printf(
" for example:\n");
3455 std::printf(
" 3 ''K width'' 1.2 0.1\n");
3456 std::printf(
" the last definition is followed by a blank line or a zero.\n");
3464 if( !strncmp(comd.c_str(),
"REL",3) ) {
3465 std::printf(
" ***>RELease <parno> [parno] ... [parno]\n");
3466 std::printf(
" If <parno> is the number of a previously variable parameter\n");
3467 std::printf(
" which has been fixed by a command: FIX <parno>, then that\n");
3468 std::printf(
" parameter will return to variable status. Otherwise a warning\n");
3469 std::printf(
" message is printed and the command is ignored.\n");
3470 std::printf(
" Note that this command operates only on parameters which were\n");
3471 std::printf(
" at one time variable and have been FIXed. It cannot make\n");
3472 std::printf(
" constant parameters variable; that must be done by redefining\n");
3473 std::printf(
" the parameter with a PARameters command.\n");
3481 if( !strncmp(comd.c_str(),
"RES",3) ) {
3482 std::printf(
" ***>REStore [code]\n");
3483 std::printf(
" If no [code] is specified, this command restores all previously\n");
3484 std::printf(
" FIXed parameters to variable status. If [code]=1, then only\n");
3485 std::printf(
" the last parameter FIXed is restored to variable status.\n");
3486 std::printf(
" If code is neither zero nor one, the command is ignored.\n");
3494 if( !strncmp(comd.c_str(),
"RET",3) ) {
3495 std::printf(
" ***>RETURN\n");
3496 std::printf(
" Signals the end of a data block, and instructs Minuit to return\n");
3497 std::printf(
" to the program which called it. The RETurn command first\n");
3498 std::printf(
" causes Minuit to CALL FCN with IFLAG=3, in order to allow FCN\n");
3499 std::printf(
" to perform any calculations associated with the final fitted\n");
3500 std::printf(
" parameter values, unless a CALL FCN 3 command has already been\n");
3501 std::printf(
" executed at the current FCN value.\n");
3509 if( !strncmp(comd.c_str(),
"SAV",3) ) {
3510 std::printf(
" ***>SAVe\n");
3511 std::printf(
" Causes the current parameter values to be saved on a file in\n");
3512 std::printf(
" such a format that they can be read in again as Minuit\n");
3513 std::printf(
" parameter definitions. If the covariance matrix exists, it is\n");
3514 std::printf(
" also output in such a format. The unit number is by default 7,\n");
3515 std::printf(
" or that specified by the user in his call to MINTIO or\n");
3516 std::printf(
" MNINIT. The user is responsible for opening the file previous\n");
3517 std::printf(
" to issuing the [SAVe] command (except where this can be done\n");
3518 std::printf(
" interactively).\n");
3526 if( !strncmp(comd.c_str(),
"SCA",3) ) {
3527 std::printf(
" ***>SCAn [parno] [numpts] [from] [to]\n");
3528 std::printf(
" Scans the value of the user function by varying parameter\n");
3529 std::printf(
" number [parno], leaving all other parameters fixed at the\n");
3530 std::printf(
" current value. If [parno] is not specified, all variable\n");
3531 std::printf(
" parameters are scanned in sequence.\n");
3532 std::printf(
" The number of points [numpts] in the scan is 40 by default,\n");
3533 std::printf(
" and cannot exceed 100. The range of the scan is by default\n");
3534 std::printf(
" 2 standard deviations on each side of the current best value,\n");
3535 std::printf(
" but can be specified as from [from] to [to].\n");
3536 std::printf(
" After each scan, if a new minimum is found, the best parameter\n");
3537 std::printf(
" values are retained as start values for future scans or\n");
3538 std::printf(
" minimizations. The curve resulting from each scan is plotted\n");
3539 std::printf(
" on the output unit in order to show the approximate behaviour\n");
3540 std::printf(
" of the function.\n");
3541 std::printf(
" This command is not intended for minimization, but is sometimes\n");
3542 std::printf(
" useful for debugging the user function or finding a\n");
3543 std::printf(
" reasonable starting point.\n");
3551 if( !strncmp(comd.c_str(),
"SEE",3) ) {
3552 std::printf(
" ***>SEEk [maxcalls] [devs]\n");
3553 std::printf(
" Causes a Monte Carlo minimization of the function, by choosing\n");
3554 std::printf(
" random values of the variable parameters, chosen uniformly\n");
3555 std::printf(
" over a hypercube centered at the current best value.\n");
3556 std::printf(
" The region size is by default 3 standard deviations on each\n");
3557 std::printf(
" side, but can be changed by specifying the value of [devs].\n");
3565 if( !strncmp(comd.c_str(),
"SET",3) ) {
3566 std::printf(
" ***>SET <option_name>\n");
3567 std::printf(
" SET BATch\n");
3568 std::printf(
" Informs Minuit that it is running in batch mode.\n");
3571 std::printf(
" SET EPSmachine <accuracy>\n");
3572 std::printf(
" Informs Minuit that the relative floating point arithmetic\n");
3573 std::printf(
" precision is <accuracy>. Minuit determines the nominal\n");
3574 std::printf(
" precision itself, but the SET EPSmachine command can be\n");
3575 std::printf(
" used to override Minuit own determination, when the user\n");
3576 std::printf(
" knows that the FCN function value is not calculated to\n");
3577 std::printf(
" the nominal machine accuracy. Typical values of <accuracy>\n");
3578 std::printf(
" are between 10**-5 and 10**-14.\n");
3581 std::printf(
" SET ERRordef <up>\n");
3582 std::printf(
" Sets the value of UP (default value= 1.), defining parameter\n");
3583 std::printf(
" errors. Minuit defines parameter errors as the change\n");
3584 std::printf(
" in parameter value required to change the function value\n");
3585 std::printf(
" by UP. Normally, for chisquared fits UP=1, and for negative\n");
3586 std::printf(
" log likelihood, UP=0.5.\n");
3589 std::printf(
" SET GRAdient [force]\n");
3590 std::printf(
" Informs Minuit that the user function is prepared to\n");
3591 std::printf(
" calculate its own first derivatives and return their values\n");
3592 std::printf(
" in the array GRAD when IFLAG=2 (see specs of FCN).\n");
3593 std::printf(
" If [force] is not specified, Minuit will calculate\n");
3594 std::printf(
" the FCN derivatives by finite differences at the current\n");
3595 std::printf(
" point and compare with the user calculation at that point,\n");
3596 std::printf(
" accepting the user values only if they agree.\n");
3597 std::printf(
" If [force]=1, Minuit does not do its own derivative\n");
3598 std::printf(
" calculation, and uses the derivatives calculated in FCN.\n");
3601 std::printf(
" SET INPut [unitno] [filename]\n");
3602 std::printf(
" Causes Minuit, in data-driven mode only, to read subsequent\n");
3603 std::printf(
" commands (or parameter definitions) from a different input\n");
3604 std::printf(
" file. If no [unitno] is specified, reading reverts to the\n");
3605 std::printf(
" previous input file, assuming that there was one.\n");
3606 std::printf(
" If [unitno] is specified, and that unit has not been opened,\n");
3607 std::printf(
" then Minuit attempts to open the file [filename]} if a\n");
3608 std::printf(
" name is specified. If running in interactive mode and\n");
3609 std::printf(
" [filename] is not specified and [unitno] is not opened,\n");
3610 std::printf(
" Minuit prompts the user to enter a file name.\n");
3611 std::printf(
" If the word REWIND is added to the command (note:no blanks\n");
3612 std::printf(
" between INPUT and REWIND), the file is rewound before\n");
3613 std::printf(
" reading. Note that this command is implemented in standard\n");
3614 std::printf(
" Fortran 77 and the results may depend on the system;\n");
3615 std::printf(
" for example, if a filename is given under VM/CMS, it must\n");
3616 std::printf(
" be preceeded by a slash.\n");
3619 std::printf(
" SET INTeractive\n");
3620 std::printf(
" Informs Minuit that it is running interactively.\n");
3623 std::printf(
" SET LIMits [parno] [lolim] [uplim]\n");
3624 std::printf(
" Allows the user to change the limits on one or all\n");
3625 std::printf(
" parameters. If no arguments are specified, all limits are\n");
3626 std::printf(
" removed from all parameters. If [parno] alone is specified,\n");
3627 std::printf(
" limits are removed from parameter [parno].\n");
3628 std::printf(
" If all arguments are specified, then parameter [parno] will\n");
3629 std::printf(
" be bounded between [lolim] and [uplim].\n");
3630 std::printf(
" Limits can be specified in either order, Minuit will take\n");
3631 std::printf(
" the smaller as [lolim] and the larger as [uplim].\n");
3632 std::printf(
" However, if [lolim] is equal to [uplim], an error condition\n");
3633 std::printf(
" results.\n");
3636 std::printf(
" SET LINesperpage\n");
3637 std::printf(
" Sets the number of lines for one page of output.\n");
3638 std::printf(
" Default value is 24 for interactive mode\n");
3641 std::printf(
" SET NOGradient");
3642 std::printf(
" The inverse of SET GRAdient, instructs Minuit not to\n");
3643 std::printf(
" use the first derivatives calculated by the user in FCN.\n");
3646 std::printf(
" SET NOWarnings\n");
3647 std::printf(
" Supresses Minuit warning messages.\n");
3650 std::printf(
" SET OUTputfile <unitno>\n");
3651 std::printf(
" Instructs Minuit to write further output to unit <unitno>.\n");
3654 std::printf(
" SET PAGethrow <integer>\n");
3655 std::printf(
" Sets the carriage control character for ``new page'' to\n");
3656 std::printf(
" <integer>. Thus the value 1 produces a new page, and 0\n");
3657 std::printf(
" produces a blank line, on some devices (see TOPofpage)\n");
3661 std::printf(
" SET PARameter <parno> <value>\n");
3662 std::printf(
" Sets the value of parameter <parno> to <value>.\n");
3663 std::printf(
" The parameter in question may be variable, fixed, or\n");
3664 std::printf(
" constant, but must be defined.\n");
3667 std::printf(
" SET PRIntout <level>\n");
3668 std::printf(
" Sets the print level, determining how much output will be\n");
3669 std::printf(
" produced. Allowed values and their meanings are displayed\n");
3670 std::printf(
" after a SHOw PRInt command, and are currently <level>=:\n");
3671 std::printf(
" [-1] no output except from SHOW commands\n");
3672 std::printf(
" [0] minimum output\n");
3673 std::printf(
" [1] default value, normal output\n");
3674 std::printf(
" [2] additional output giving intermediate results.\n");
3675 std::printf(
" [3] maximum output, showing progress of minimizations.\n");
3676 std::printf(
" Note: See also the SET WARnings command.\n");
3679 std::printf(
" SET RANdomgenerator <seed>\n");
3680 std::printf(
" Sets the seed of the random number generator used in SEEk.\n");
3681 std::printf(
" This can be any integer between 10000 and 900000000, for\n");
3682 std::printf(
" example one which was output from a SHOw RANdom command of\n");
3683 std::printf(
" a previous run.\n");
3686 std::printf(
" SET STRategy <level>\n");
3687 std::printf(
" Sets the strategy to be used in calculating first and second\n");
3688 std::printf(
" derivatives and in certain minimization methods.\n");
3689 std::printf(
" In general, low values of <level> mean fewer function calls\n");
3690 std::printf(
" and high values mean more reliable minimization.\n");
3691 std::printf(
" Currently allowed values are 0, 1 (default), and 2.\n");
3694 std::printf(
" SET TITle\n");
3695 std::printf(
" Informs Minuit that the next input line is to be considered\n");
3696 std::printf(
" the (new) title for this task or sub-task. This is for\n");
3697 std::printf(
" the convenience of the user in reading his output.\n");
3700 std::printf(
" SET WARnings\n");
3701 std::printf(
" Instructs Minuit to output warning messages when suspicious\n");
3702 std::printf(
" conditions arise which may indicate unreliable results.\n");
3703 std::printf(
" This is the default.\n");
3706 std::printf(
" SET WIDthpage\n");
3707 std::printf(
" Informs Minuit of the output page width.\n");
3708 std::printf(
" Default values are 80 for interactive jobs\n");
3716 if( !strncmp(comd.c_str(),
"SHO",3) ) {
3717 std::printf(
" ***>SHOw <option_name>\n");
3718 std::printf(
" All SET XXXX commands have a corresponding SHOw XXXX command.\n");
3719 std::printf(
" In addition, the SHOw commands listed starting here have no\n");
3720 std::printf(
" corresponding SET command for obvious reasons.\n");
3723 std::printf(
" SHOw CORrelations\n");
3724 std::printf(
" Calculates and prints the parameter correlations from the\n");
3725 std::printf(
" error matrix.\n");
3728 std::printf(
" SHOw COVariance\n");
3729 std::printf(
" Prints the (external) covariance (error) matrix.\n");
3732 std::printf(
" SHOw EIGenvalues\n");
3733 std::printf(
" Calculates and prints the eigenvalues of the covariance\n");
3734 std::printf(
" matrix.");
3737 std::printf(
" SHOw FCNvalue\n");
3738 std::printf(
" Prints the current value of FCN.\n");
3746 if( !strncmp(comd.c_str(),
"SIM",3) ) {
3747 std::printf(
" ***>SIMplex [maxcalls] [tolerance]\n");
3748 std::printf(
" Performs a function minimization using the simplex method of\n");
3749 std::printf(
" Nelder and Mead. Minimization terminates either when the\n");
3750 std::printf(
" function has been called (approximately) [maxcalls] times,\n");
3751 std::printf(
" or when the estimated vertical distance to minimum (EDM) is\n");
3752 std::printf(
" less than [tolerance].\n");
3753 std::printf(
" The default value of [tolerance] is 0.1*UP(see SET ERRordef).\n");
3761 if( !strncmp(comd.c_str(),
"STA",3) ) {
3762 std::printf(
" ***>STAndard\n");
3770 if( !strncmp(comd.c_str(),
"STO",3) ) {
3771 std::printf(
" ***>STOP\n");
3772 std::printf(
" Same as EXIT.\n");
3780 if( !strncmp(comd.c_str(),
"TOP",3) ) {
3781 std::printf(
" ***>TOPofpage\n");
3782 std::printf(
" Causes Minuit to write the character specified in a\n");
3783 std::printf(
" SET PAGethrow command (default = 1) to column 1 of the output\n");
3784 std::printf(
" file, which may or may not position your output medium to\n");
3785 std::printf(
" the top of a page depending on the device and system.\n");
3789 std::printf(
" Unknown MINUIT command. Type HELP for list of commands.\n");
3808 Double_urt dmin_, dxdi, elem, wint, tlrg2, d, dlast, ztemp, g2bfor;
3809 Double_urt df, aimsag, fs1, tlrstp, fs2, stpinm, g2i, sag=0, xtf, xti, xtj;
3810 Int_urt icyc, ncyc, ndex, idrv, iext, npar2, i, j, ifail, npard, nparx, id, multpy;
3813 ldebug =
fIdbg[3] >= 1;
3830 if (
fISW[4] >= 2 || ldebug) {
3831 std::printf(
" START COVARIANCE MATRIX CALCULATION.\n");
3844 ostringstream warning;
3845 warning <<
"function value differs from AMIN by " << df;
3846 mnwarn(
"D",
"MNHESS", warning.str().c_str());
3850 std::printf(
" PAR D GSTEP D G2 GRD SAG \n");
3860 for (i = 1; i <= npar2; ++i) {
fVhmat[i-1] = 0; }
3864 for (
id = 1;
id <= npard; ++id) {
3865 i =
id +
fNpar - npard;
3867 if (
fG2[i-1] == 0) {
3868 ostringstream warning;
3869 warning <<
"Second derivative enters zero, param " << iext;
3870 mnwarn(
"W",
"HESSE", warning.str().c_str());
3873 if (m_userParameterFlag[iext] > 1) {
3878 fG2[i-1] =
fUp / (wint*wint);
3886 for (icyc = 1; icyc <= ncyc; ++icyc) {
3888 for (multpy = 1; multpy <= 5; ++multpy) {
3900 sag = (fs1 + fs2 - fAmin*2)*.5;
3901 if (sag != 0)
goto L30;
3903 if (d >= .5)
goto L26;
3905 if (d > .5) d = .51;
3911 { ostringstream warning3;
3912 warning3 <<
"Second derivative zero for parameter " << iext;
3913 mnwarn(
"W",
"HESSE", warning3.str().c_str()); }
3918 fG2[i-1] = sag*2 / (d*d);
3919 fGrd[i-1] = (fs1 - fs2) / (d*2);
3921 std::printf(
"%4d%2d%12.5g%12.5g%12.5g%12.5g\n",i,idrv,
fGstep[i-1],
fG2[i-1],
fGrd[i-1],sag);
3932 if (d < dmin_) d = dmin_;
3934 if (URMath::Abs((d - dlast) / d) < tlrstp ||
3935 URMath::Abs((
fG2[i-1] - g2bfor) /
fG2[i-1]) < tlrg2) {
3944 ostringstream warning;
3945 warning <<
"Second Deriv. SAG,AIM= " << iext <<
" " << sag <<
" " << aimsag;
3946 mnwarn(
"D",
"MNHESS", warning.str().c_str());
3948 ndex = i*(i + 1) / 2;
3959 if (
fNpar == 1)
goto L214;
3960 for (i = 1; i <=
fNpar; ++i) {
3961 for (j = 1; j <= i-1; ++j) {
3965 fX[j-1] = xtj + fDirin[j-1];
3972 fDirin[i-1]*fDirin[j-1]);
3973 ndex = i*(i-1) / 2 + j;
3981 for (i = 1; i <=
fNpar; ++i) {
3982 for (j = 1; j <= i; ++j) {
3983 ndex = i*(i-1) / 2 + j;
3994 mnwarn(
"W",
"HESSE",
"Matrix inversion fails.");
4000 for (i = 1; i <=
fNpar; ++i) {
4003 for (j = 1; j <= i-1; ++j) {
4015 std::printf(
" COVARIANCE MATRIX CALCULATED SUCCESSFULLY\n");
4024 std::printf(
" MNHESS FAILS AND WILL RETURN DIAGONAL MATRIX. \n");
4026 for (i = 1; i <=
fNpar; ++i) {
4028 for (j = 1; j <= i-1; ++j) {
4034 if (g2i <= 0) g2i = 1;
4035 fVhmat[ndex-1] = 2 / g2i;
4051 Double_urt dmin_, d, dfmin, dgmin=0, change, chgold, grdold=0, epspri;
4052 Double_urt fs1, optstp, fs2, grdnew=0, sag, xtf;
4053 Int_urt icyc, ncyc=0, idrv, i, nparx;
4056 ldebug =
fIdbg[5] >= 1;
4064 for (i = 1; i <=
fNpar; ++i) {
4070 if (d > optstp) d = optstp;
4071 if (d < dmin_) d = dmin_;
4074 for (icyc = 1; icyc <= ncyc; ++icyc) {
4085 sag = (fs1 + fs2 -
fAmin*2)*.5;
4087 grdnew = (fs1 - fs2) / (d*2);
4090 std::printf(
"%4d%2d%12.5g%12.5g%12.5g%12.5g%12.5g\n",i,idrv,
fGstep[i-1],d,
fG2[i-1],grdnew,sag);
4092 if (grdnew == 0)
goto L60;
4094 if (change > chgold && icyc > 1)
goto L60;
4100 if (change < .05)
goto L60;
4101 if (URMath::Abs(grdold - grdnew) < dgmin)
goto L60;
4103 mnwarn(
"D",
"MNHES1",
"Step size too small for 1st drv.");
4109 { std::ostringstream warning;
4110 warning <<
"Too many iterations on D1. " << grdold <<
" " << grdnew;
4111 mnwarn(
"D",
"MNHES1", warning.str().c_str()); }
4137 Double_urt pb, ep, wg, xi, sigsav, reg, sig2;
4138 Int_urt npfn, ndex, loop=0, i, j, ifail, iseed;
4139 Int_urt jhold, nloop, nparx, nparp1, jh, jl, iswtr;
4141 if (
fNpar <= 0)
return;
4148 if (nloop <= 0) nloop =
fNpar + 4;
4155 for (i = 1; i <=
fNpar; ++i) {
4158 for (j = 1; j <= i; ++j) {
4159 ndex = i*(i-1) / 2 + j;
4165 if (ifail >= 1)
goto L280;
4167 for (i = 1; i <=
fNpar; ++i) {
4169 for (j = 1; j <= i; ++j) {
4177 for (i = 1; i <=
fNpar; ++i) {
4185 std::printf(
"START ATTEMPT NO.%2d TO FIND NEW MINIMUM\n",loop);
4195 for (i = 1; i <=
fNpar; ++i) {
4198 fX[i-1] = xi -
fDirin[i-1]*(rnum - .5);
4204 }
else if (
fIMPRy[i-1] > amax) {
4217 if (
fAmin < 0)
goto L95;
4218 if (
fISW[1] <= 2)
goto L280;
4220 if (sig2 < ep &&
fEDM < ep)
goto L100;
4224 for (i = 1; i <=
fNpar; ++i) {
4232 if (ystar >=
fAmin)
goto L70;
4234 for (i = 1; i <=
fNpar; ++i) {
4239 if (ystst <
fIMPRy[jl-1])
goto L67;
4247 if (ystar >=
fIMPRy[jh-1])
goto L73;
4250 if (jhold != jh)
goto L50;
4253 for (i = 1; i <=
fNpar; ++i) {
4258 if (ystst >
fIMPRy[jh-1])
goto L30;
4260 if (ystst <
fAmin)
goto L67;
4266 std::printf(
" AN IMPROVEMENT ON THE PREVIOUS MINIMUM HAS BEEN FOUND\n");
4274 for (i = 1; i <=
fNpar; ++i) {
4284 for (i = 1; i <=
fNpar; ++i) {
4302 std::printf(
" IMPROVE HAS FOUND A TRULY NEW MINIMUM\n");
4303 std::printf(
" *************************************\n");
4309 std::printf(
" COVARIANCE MATRIX WAS NOT POSITIVE-DEFINITE\n");
4315 for (i = 1; i <=
fNpar; ++i) {
4324 std::printf(
" IMPROVE HAS RETURNED TO REGION OF ORIGINAL MINIMUM\n");
4328 if (
fISW[1] < 2)
goto L380;
4329 if (loop < nloop &&
fISW[0] < 1)
goto L20;
4346 for (j = 1; j <=
fNpar; ++j) {
4349 if (m_userParameterFlag[i] == 1) {
4351 m_userParameterValue[i] = pint[j-1];
4369 Double_urt piby2, epsp1, epsbak, epstry, distnn;
4387 fCovmes[0] =
"NO ERROR MATRIX ";
4388 fCovmes[1] =
"ERR MATRIX APPROXIMATE";
4389 fCovmes[2] =
"ERR MATRIX NOT POS-DEF";
4390 fCovmes[3] =
"ERROR MATRIX ACCURATE ";
4408 for (idb = 0; idb <= 10; ++idb) {
fIdbg[idb] = 0; }
4428 for (i = 1; i <= 100; ++i) {
4432 if (epsbak < epstry)
goto L35;
4436 std::printf(
" MNINIT UNABLE TO DETERMINE ARITHMETIC PRECISION. WILL ASSUME:%g\n",
fEpsmac);
4461 Int_urt kint, i2, newcod, ifx=0, inu;
4467 if (i2 >
fMaxext || i2 < 0)
goto L900;
4468 if (i2 > 0)
goto L30;
4472 for (inu = 1; inu <=
fNu; ++inu) {
4474 if (m_userParameterFlag[inu] <= 0)
continue;
4476 if (m_userParameterFlag[inu] == 1 && newcod == 1)
continue;
4478 kint = m_userParameterIdToInternalId[inu];
4482 std::printf(
" LIMITS NOT CHANGED FOR FIXED PARAMETER:%4d\n",inu);
4489 std::printf(
" LIMITS REMOVED FROM PARAMETER :%3d\n",inu);
4493 snew =
fGstep[kint-1]*dxdi;
4496 m_userParameterFlag[inu] = 1;
4502 std::printf(
" PARAMETER %3d LIMITS SET TO %15.5g%15.5g\n",inu,
fAlim[inu-1],
fBlim[inu-1]);
4505 m_userParameterFlag[inu] = 4;
4514 if (m_userParameterFlag[i2] <= 0) {
4515 std::printf(
" PARAMETER %3d IS NOT VARIABLE.\n", i2);
4519 kint = m_userParameterIdToInternalId[i2];
4522 std::printf(
" REQUEST TO CHANGE LIMITS ON FIXED PARAMETER:%3d\n",i2);
4523 for (ifx = 1; ifx <=
fNpfix; ++ifx) {
4524 if (i2 ==
fIpfix[ifx-1])
goto L92;
4526 std::printf(
" MINUIT BUG IN MNLIMS. SEE F. JAMES\n");
4533 if (m_userParameterFlag[i2] != 1) {
4535 std::printf(
" LIMITS REMOVED FROM PARAMETER %2d\n",i2);
4544 fGrd[kint-1] *= dxdi;
4547 m_userParameterFlag[i2] = 1;
4549 std::printf(
" NO LIMITS SPECIFIED. PARAMETER %3d IS ALREADY UNLIMITED. NO CHANGE.\n",i2);
4557 m_userParameterFlag[i2] = 4;
4559 std::printf(
" PARAMETER %3d LIMITS SET TO %15.5g%15.5g\n",i2,
fAlim[i2-1],
fBlim[i2-1]);
4562 if (kint <= 0)
fGsteps[ifx-1] = -.1;
4563 else fGstep[kint-1] = -.1;
4566 if (
fCstatu !=
"NO CHANGE ") {
4591 Double_urt xpq[12], ypq[12], slam, sdev, coeff[3], denom, flast;
4592 Double_urt fvals[3], xvals[3], f1, fvmin, xvmin, ratio, f2, f3, fvmax;
4593 Double_urt toler8, toler9, overal, undral, slamin, slamax, slopem;
4594 Int_urt i, nparx=0, nvmax=0, nxypt, kk, ipt;
4601 l65 = 0; l70 = 0; l80 = 0;
4602 ldebug =
fIdbg[1] >= 1;
4612 std::printf(
" MNLINE start point not consistent, F values, parameters=\n");
4613 for (kk = 1; kk <=
fNpar; ++kk) {
4614 std::printf(
" %14.5e\n",
fX[kk-1]);
4627 for (i = 1; i <=
fNpar; ++i) {
4628 if (step[i-1] != 0) {
4630 if (slamin == 0) slamin = ratio;
4631 if (ratio < slamin) slamin = ratio;
4633 fX[i-1] = start[i-1] + step[i-1];
4635 if (slamin == 0) slamin =
fEpsmac;
4643 chpq[nxypt-1] =
charal[nxypt-1];
4658 denom = (flast - fstart - slope*slam)*2 / (slam*slam);
4660 if (denom != 0) slam = -slope / denom;
4661 if (slam < 0) slam = slamax;
4662 if (slam > slamax) slam = slamax;
4663 if (slam < toler8) slam = toler8;
4664 if (slam < slamin) {
4668 if (
URMath::Abs(slam - 1) < toler8 && f1 < fstart) {
4672 if (
URMath::Abs(slam - 1) < toler8) slam = toler8 + 1;
4677 for (i = 1; i <=
fNpar; ++i) {
fX[i-1] = start[i-1] + slam*step[i-1]; }
4683 chpq[nxypt-1] =
charal[nxypt-1];
4684 xpq[nxypt-1] = slam;
4690 if (fstart == fvmin) {
4692 toler8 = toler*slam;
4693 overal = slam - toler8;
4696 }
while (fstart == fvmin);
4698 if (!l65 && !l70 && !l80) {
4702 xvals[1] = xpq[nxypt-2];
4703 fvals[1] = ypq[nxypt-2];
4704 xvals[2] = xpq[nxypt-1];
4705 fvals[2] = ypq[nxypt-1];
4709 mnpfit(xvals, fvals, 3, coeff, sdev);
4710 if (coeff[2] <= 0) {
4711 slopem = coeff[2]*2*xvmin + coeff[1];
4712 if (slopem <= 0) slam = xvmin + slamax;
4713 else slam = xvmin - slamax;
4715 slam = -coeff[1] / (coeff[2]*2);
4716 if (slam > xvmin + slamax) slam = xvmin + slamax;
4717 if (slam < xvmin - slamax) slam = xvmin - slamax;
4719 if (slam > 0) {
if (slam > overal) slam = overal; }
4720 else {
if (slam < undral) slam = undral; }
4725 for (ipt = 1; ipt <= 3; ++ipt) {
4737 for (i = 1; i <=
fNpar; ++i) {
fX[i-1] = start[i-1] + slam*step[i-1]; }
4742 chpq[nxypt-1] =
charal[nxypt-1];
4743 xpq[nxypt-1] = slam;
4748 if (fvals[1] > fvmax) {
4752 if (fvals[2] > fvmax) {
4762 if (slam > xvmin) overal =
URMath::Min(overal,slam - toler8);
4763 if (slam < xvmin) undral =
URMath::Max(undral,slam + toler8);
4764 slam = (slam + xvmin)*.5;
4766 }
while (f3 >= fvmax);
4769 if (l65 || l70)
break;
4771 xvals[nvmax-1] = slam;
4772 fvals[nvmax-1] = f3;
4777 if (slam > xvmin) overal =
URMath::Min(overal,slam - toler8);
4778 if (slam < xvmin) undral =
URMath::Max(undral,slam + toler8);
4780 }
while (nxypt < 12);
4786 cmess =
" LINE SEARCH HAS EXHAUSTED THE LIMIT OF FUNCTION CALLS ";
4788 std::printf(
" MNLINE DEBUG: steps=\n");
4789 for (kk = 1; kk <=
fNpar; ++kk) {
4790 std::printf(
" %12.4g\n",step[kk-1]);
4795 if (l70) cmess =
" LINE SEARCH HAS ATTAINED TOLERANCE ";
4796 if (l80) cmess =
" STEP SIZE AT ARITHMETICALLY ALLOWED MINIMUM";
4799 for (i = 1; i <=
fNpar; ++i) {
4800 fDirin[i-1] = step[i-1]*xvmin;
4801 fX[i-1] = start[i-1] +
fDirin[i-1];
4805 mnwarn(
"D",
"MNLINE",
" LINE MINIMUM IN BACKWARDS DIRECTION");
4807 if (fvmin == fstart) {
4808 mnwarn(
"D",
"MNLINE",
" LINE SEARCH FINDS NO IMPROVEMENT ");
4811 std::printf(
" AFTER%3d POINTS,%s\n",nxypt,cmess.c_str());
4826 Int_urt ndex, i, j, m, n, ncoef, nparm, id, it, ix;
4827 Int_urt nsofar, ndi, ndj, iso, isw2, isw5;
4832 std::printf(
"%s\n",(
fCovmes[isw2]).c_str());
4836 std::printf(
" MNMATU: NPAR=0\n");
4845 std::printf(
"%s\n",(
fCovmes[isw2]).c_str());
4850 if (
fNpar <= 1)
return;
4856 std::printf(
" PARAMETER CORRELATION COEFFICIENTS \n");
4857 ctemp =
" NO. GLOBAL";
4858 for (
id = 1;
id <= nparm; ++id) {
4859 ostringstream localTemp;
4860 localTemp << setw(7) << right <<
fNexofi[
id-1];
4861 ctemp += localTemp.str();
4863 std::printf(
"%s\n",ctemp.c_str());
4864 for (i = 1; i <=
fNpar; ++i) {
4866 ndi = i*(i + 1) / 2;
4867 for (j = 1; j <=
fNpar; ++j) {
4870 ndex = m*(m-1) / 2 + n;
4871 ndj = j*(j + 1) / 2;
4875 ostringstream octemp;
4876 octemp << setw(7) << right << ix <<
": " << fixed << setw(9) << setprecision(3) <<
fGlobcc[i-1];
4877 ctemp = octemp.str();
4878 for (it = 1; it <= nparm; ++it) {
4879 ostringstream localTemp;
4880 localTemp << fixed << setw(7) << setprecision(3) <<
fMATUvline[it-1];
4881 ctemp += localTemp.str();
4883 std::printf(
"%s\n",ctemp.c_str());
4884 if (i <= nparm)
continue;
4885 for (iso = 1; iso <= 10; ++iso) {
4889 for (it = nsofar + 1; it <= nparm; ++it) {
4890 ostringstream localTemp;
4891 localTemp << fixed << setw(7) << setprecision(3) <<
fMATUvline[it-1];
4892 ctemp = ctemp + localTemp.str();
4894 std::printf(
"%s\n",ctemp.c_str());
4895 if (i <= nparm)
break;
4899 std::printf(
" %s\n",(
fCovmes[isw2]).c_str());
4914 Double_urt gdel, gami, vlen, dsum, gssq, vsum, d;
4917 Int_urt npfn, ndex, iext, i, j, m, n, npsdf, nparx;
4918 Int_urt iswtr, lined2, kk, nfcnmg, nrstrt,iter;
4922 if (
fNpar <= 0)
return;
4937 rhotol =
fApsi*.001;
4939 std::printf(
" START MIGRAD MINIMIZATION. STRATEGY%2d. CONVERGENCE WHEN EDM .LT.%9.2e\n",
fIstrat,rhotol);
4954 if (
fISW[1] >= 1)
goto L10;
4963 if (
fISW[1] >= 1)
goto L10;
4965 for (i = 1; i <=
fNpar; ++i) {
4972 for (i = 1; i <=
fNpar; ++i) {
4973 if (
fG2[i-1] > 0)
continue;
4979 mnwarn(
"D",
"MNMIGR",
"Negative G2 line search");
4982 std::printf(
" Negative G2 line search, param %3d %13.3g%13.3g\n",iext,fs,
fAmin);
4988 for (i = 1; i <=
fNpar; ++i) {
4990 for (j = 1; j <= i-1; ++j) {
4995 if (
fG2[i-1] <= 0)
fG2[i-1] = 1;
5000 std::printf(
" DEBUG MNMIGR, STARTING MATRIX DIAGONAL, VHMAT=\n");
5001 for (kk = 1; kk <=
Int_urt(vlen); ++kk) {
5002 std::printf(
" %10.2g\n",
fVhmat[kk-1]);
5015 for (i = 1; i <=
fNpar; ++i) {
5019 for (j = 1; j <= i-1; ++j) {
5028 mnwarn(
"W",
"MIGRAD",
"STARTING MATRIX NOT POS-DEFINITE.");
5038 if (iswtr >= 2)
mnmatu(0);
5044 for (i = 1; i <=
fNpar; ++i) {
5047 for (j = 1; j <=
fNpar; ++j) {
5050 ndex = m*(m-1) / 2 + n;
5051 ri +=
fVhmat[ndex-1]*fMIGRgs[j-1];
5057 mnwarn(
"D",
"MIGRAD",
" FIRST DERIVATIVES OF FCN ARE ALL ZERO");
5062 mnwarn(
"D",
"MIGRAD",
" NEWTON STEP NOT DESCENT.");
5063 if (npsdf == 1)
goto L1;
5070 if (
fAmin == fs)
goto L200;
5088 for (i = 1; i <=
fNpar; ++i) {
5091 for (j = 1; j <=
fNpar; ++j) {
5094 ndex = m*(m-1) / 2 + n;
5100 gdgssq += gami*gami;
5102 delgam +=
fDirin[i-1]*gami;
5107 if (
fEDM < 0 || gvg <= 0) {
5108 mnwarn(
"D",
"MIGRAD",
"NOT POS-DEF. EDM OR GVG NEGATIVE.");
5110 if (npsdf == 1)
goto L230;
5117 if (iswtr >= 3 || ( iswtr == 2 && iter % 10 == 1 ) ) {
5122 mnwarn(
"D",
"MIGRAD",
"NO CHANGE IN FIRST DERIVATIVES OVER LAST STEP");
5125 mnwarn(
"D",
"MIGRAD",
"FIRST DERIVATIVES INCREASING ALONG SEARCH LINE");
5130 std::printf(
" VHMAT 1 =\n");
5131 for (kk = 1; kk <= 10; ++kk) {
5132 std::printf(
" %10.2g\n",
fVhmat[kk-1]);
5137 for (i = 1; i <=
fNpar; ++i) {
5138 for (j = 1; j <= i; ++j) {
5139 if(delgam == 0 || gvg == 0) d = 0;
5142 ndex = i*(i-1) / 2 + j;
5149 if (iswtr >= 3 || ldebug) {
5150 std::printf(
" RELATIVE CHANGE IN COV. MATRIX=%5.1f per cent\n",
fDcovar*100);
5153 std::printf(
" VHMAT 2 =\n");
5154 for (kk = 1; kk <= 10; ++kk) {
5155 std::printf(
" %10.3g\n",
fVhmat[kk-1]);
5158 if (delgam <= gvg)
goto L135;
5159 for (i = 1; i <=
fNpar; ++i) {
5162 for (i = 1; i <=
fNpar; ++i) {
5163 for (j = 1; j <= i; ++j) {
5164 ndex = i*(i-1) / 2 + j;
5170 if (
fEDM < rhotol*.1)
goto L300;
5172 for (i = 1; i <=
fNpar; ++i) {
5186 std::printf(
" CALL LIMIT EXCEEDED IN MIGRAD.\n");
5193 std::printf(
" MIGRAD FAILS TO FIND IMPROVEMENT\n");
5196 if (
fEDM < rhotol)
goto L300;
5199 std::printf(
" MACHINE ACCURACY LIMITS FURTHER IMPROVEMENT.\n");
5205 std::printf(
" MIGRAD FAILS WITH STRATEGY=0. WILL TRY WITH STRATEGY=1.\n");
5213 std::printf(
" MIGRAD TERMINATED WITHOUT CONVERGENCE.\n");
5221 std::printf(
" MIGRAD MINIMIZATION HAS CONVERGED.\n");
5226 std::printf(
" MIGRAD WILL VERIFY CONVERGENCE AND ERROR MATRIX.\n");
5231 if (
fEDM > rhotol)
goto L10;
5242 if (iswtr >= 0)
mnprin(3, fAmin);
5243 if (iswtr >= 1)
mnmatu(1);
5258 Int_urt nbad, ilax, ilax2, ngood, nfcnmi, iin, knt;
5260 if (
fNpar <= 0)
goto L700;
5265 for (knt = 1; knt <=
fNpar; ++knt) {
5269 if (knt >= 7)
break;
5271 if (ilax == 0)
break;
5272 if (ilax > 0 && ilax <=
fNu) {
5274 if (m_userParameterIdToInternalId[ilax] > 0)
goto L565;
5276 std::printf(
" PARAMETER NUMBER %3d NOT VARIABLE. IGNORED.\n",ilax);
5282 mnmnot(ilax, ilax2, val2pl, val2mi);
5286 iin = m_userParameterIdToInternalId[ilax];
5287 if (
fErp[iin-1] > 0) ++ngood;
5289 if (
fErn[iin-1] < 0) ++ngood;
5297 if (ngood == 0 && nbad == 0)
goto L700;
5298 if (ngood > 0 && nbad == 0)
fCstatu =
"SUCCESSFUL";
5299 if (ngood == 0 && nbad > 0)
fCstatu =
"FAILURE ";
5300 if (ngood > 0 && nbad > 0)
fCstatu =
"PROBLEMS ";
5310 std::printf(
" NEW MINIMUM FOUND. GO BACK TO MINIMIZATION STEP.\n");
5311 std::printf(
" =================================================\n");
5312 std::printf(
" V\n");
5313 std::printf(
" V\n");
5314 std::printf(
" V\n");
5315 std::printf(
" VVVVVVV\n");
5316 std::printf(
" VVVVV\n");
5317 std::printf(
" VVV\n");
5318 std::printf(
" V\n");
5322 std::printf(
" THERE ARE NO MINOS ERRORS TO CALCULATE.\n");
5340 Double_urt abest, xunit, dc, ut, sigsav, du1;
5342 Int_urt marc, isig, mpar, ndex, imax, indx, ierr, i, j;
5343 Int_urt iercr, it, istrav, nfmxin, nlimit, isw2, isw4;
5357 for (i = 1; i <= mpar; ++i) {
fXt[i-1] =
fX[i-1]; }
5358 i__1 = mpar*(mpar + 1) / 2;
5360 for (i = 1; i <= mpar; ++i) {
5365 it = m_userParameterIdToInternalId[ilax];
5371 ut = m_userParameterValue[ilax];
5372 if (m_userParameterFlag[ilax] == 1) {
5374 fBlim[ilax-1] = ut + fMNOTw[it-1]*100;
5376 ndex = it*(it + 1) / 2;
5379 for (i = 1; i <= mpar; ++i) {
5380 if (i == it)
continue;
5389 std::printf(
" MINUIT ERROR. CANNOT FIX PARAMETER%4d INTERNAL%3d\n",ilax,it);
5395 for (isig = 1; isig <= 2; ++isig) {
5409 std::printf(
" DETERMINATION OF %sTIVE MINOS ERROR FOR PARAMETER%d (%s)\n" 5411 ,(m_userParameterName[ilax]).c_str());
5414 mnwarn(
"D",
"MINOS",
"NO COVARIANCE MATRIX.");
5416 nlimit =
fNfcn + nfmxin;
5422 m_userParameterValue[ilax] = ut + sig*du1;
5423 m_userParameterValue[ilax] =
URMath::Min(m_userParameterValue[ilax],
fBlim[ilax-1]);
5424 m_userParameterValue[ilax] =
URMath::Max(m_userParameterValue[ilax],
fAlim[ilax-1]);
5426 delu = m_userParameterValue[ilax] - ut;
5430 fac = delu /
fMNOTw[it-1];
5431 for (i = 1; i <=
fNpar; ++i) {
5436 std::printf(
" PARAMETER%4d SET TO%11.3e + %10.3e = %12.3e\n",ilax,ut,delu,m_userParameterValue[ilax]);
5442 fXmidcr = m_userParameterValue[ilax];
5448 if (abest -
fAmin >
fUp*.01)
goto L650;
5449 if (iercr == 1)
goto L440;
5450 if (iercr == 2)
goto L450;
5451 if (iercr == 3)
goto L460;
5458 std::printf(
" THE %4sTIVE MINOS ERROR OF PARAMETER%3d %10s, IS %12.4e\n" 5460 ,(m_userParameterName[ilax]).c_str(),eros);
5469 std::printf(
" THE %4sTIVE MINOS ERROR OF PARAMETER%3d, %s EXCEEDS ITS LIMIT.\n" 5471 ,(m_userParameterName[ilax]).c_str());
5477 std::printf(
" THE %4sTIVE MINOS ERROR%4d REQUIRES MORE THAN%5d FUNCTION CALLS.\n" 5478 ,csig.c_str(),ilax,nfmxin);
5484 std::printf(
" %4sTIVE MINOS ERROR NOT CALCULATED FOR PARAMETER%d\n" 5485 ,csig.c_str(),ilax);
5491 std::printf(
" **************************************************************************\n");
5496 if (ilax2 > 0 && ilax2 <=
fNu) val2mi = m_userParameterValue[ilax2];
5500 if (ilax2 > 0 && ilax2 <=
fNu) val2pl = m_userParameterValue[ilax2];
5507 i__1 = mpar*(mpar + 1) / 2;
5509 for (i = 1; i <= mpar; ++i) {
5528 sav = m_userParameterValue[ilax];
5532 m_userParameterValue[ilax] = sav;
5560 Double_urt vplu, a_small, gsmin, pinti, vminu, danger, sav, sav2;
5561 Int_urt ierr, kint, in, ix, ktofix, lastin, kinfix, nvl;
5562 string cnamk, chbufi;
5571 *m_logStream <<
" MINUIT USER ERROR. PARAMETER NUMBER IS " << setw(3) << k
5572 <<
" ALLOWED RANGE IS ONE TO " << setw(4) << fMaxext <<
'\n';
5578 if (m_userParameterFlag[k] < 0)
goto L50;
5581 for (ix = 1; ix <=
fNpfix; ++ix) {
5582 if (
fIpfix[ix-1] == k) ktofix = k;
5585 mnwarn(
"W",
"PARAM DEF",
"REDEFINING A FIXED PARAMETER.");
5588 *m_logStream <<
" CANNOT RELEASE. MAX NPAR EXCEEDED.\n";
5595 if (m_userParameterIdToInternalId[k] > 0) kint =
fNpar - 1;
5602 *m_logStream <<
" PARAMETER DEFINITIONS:\n";
5603 *m_logStream <<
" NO. NAME VALUE STEP SIZE LIMITS\n";
5606 if (wk > 0)
goto L122;
5610 *m_logStream << setw(6) << k <<
" " << setw(20) << cnamk << setw(13) << setprecision(5) << uk <<
" constant\n";
5615 if (a == 0 && b == 0) {
5620 *m_logStream << setw(6) << k <<
" " << setw(20) << cnamk
5621 << setw(13) << setprecision(5) << uk
5622 << setw(13) << setprecision(5) << wk <<
" no limits\n";
5630 *m_logStream << setw(6) << k <<
" " << setw(20) << cnamk
5631 << setw(13) << setprecision(5) << uk
5632 << setw(13) << setprecision(5) << wk
5633 << setw(13) << setprecision(5) << a
5634 << setw(13) << setprecision(5) << b <<
'\n';
5641 *m_logStream <<
" MINUIT USER ERROR. TOO MANY VARIABLE PARAMETERS.\n";
5644 if (nvl == 1)
goto L200;
5646 std::printf(
" USER ERROR IN MINUIT PARAMETER\n");
5647 std::printf(
" DEFINITION\n");
5648 std::printf(
" UPPER AND LOWER LIMITS EQUAL.\n");
5649 *m_logStream <<
" USER ERROR IN MINUIT PARAMETER\n" <<
" DEFINITION\n" <<
" UPPER AND LOWER LIMITS EQUAL.\n";
5656 mnwarn(
"W",
"PARAM DEF",
"PARAMETER LIMITS WERE REVERSED.");
5660 ostringstream warning;
5661 warning <<
"LIMITS ON PARAM " << k <<
" TOO FAR APART.";
5662 mnwarn(
"W",
"PARAM DEF", warning.str().c_str());
5665 danger = (b - uk)*(uk - a);
5667 mnwarn(
"W",
"PARAM DEF",
"STARTING VALUE OUTSIDE LIMITS.");
5670 mnwarn(
"W",
"PARAM DEF",
"STARTING VALUE IS AT LIMIT.");
5681 m_userParameterName[k] = cnamk;
5682 m_userParameterValue[k] = uk;
5686 m_userParameterFlag[k] = nvl;
5692 for (ix = 1; ix <= k-1; ++ix) {
if (m_userParameterIdToInternalId[ix] > 0) ++lastin; }
5694 if (kint ==
fNpar)
goto L280;
5697 for (in =
fNpar; in >= lastin + 1; --in) {
5700 m_userParameterIdToInternalId[ix] = in + 1;
5710 for (in = lastin + 1; in <= kint; ++in) {
5713 m_userParameterIdToInternalId[ix] = in;
5725 m_userParameterIdToInternalId[ix] = 0;
5732 m_userParameterIdToInternalId[ix] = in;
5735 sav = m_userParameterValue[ix];
5738 fXt[in-1] =
fX[in-1];
5743 vplu = pinti -
fX[in-1];
5747 vminu = pinti - fX[in-1];
5759 if (m_userParameterFlag[k] > 1) {
5766 kinfix = m_userParameterIdToInternalId[ktofix];
5767 if (kinfix > 0)
mnfixp(kinfix-1, ierr);
5768 if (ierr > 0)
goto L800;
5772 *m_logStream << endl;
5776 *m_logStream << endl;
5797 Int_urt k, llist, ibegin, lenbuf, istart, lnc, icy;
5798 string cnamk, comand, celmnt, ctemp;
5801 lenbuf = strlen(crdbuf.c_str());
5803 kapo1 = strspn(crdbuf.c_str(),
"'");
5804 if (kapo1 == 0)
goto L150;
5805 kapo2 = strspn(crdbuf.c_str() + kapo1,
"'");
5806 if (kapo2 == 0)
goto L150;
5810 for (istart = 1; istart <= kapo1-1; ++istart) {
5811 if (crdbuf.substr(istart-1,1) !=
" ")
goto L120;
5816 celmnt = crdbuf.substr(istart-1, kapo1-istart);
5817 scanf(celmnt.c_str(),fk);
5819 if (k <= 0)
goto L210;
5820 cnamk =
"PARAM " + celmnt;
5821 if (kapo2 - kapo1 > 1) {
5822 cnamk = crdbuf.substr(kapo1, kapo2-1-kapo1);
5825 for (icy = kapo2 + 1; icy <= lenbuf; ++icy) {
5826 if (crdbuf.substr(icy-1,1) ==
",")
goto L139;
5827 if (crdbuf.substr(icy-1,1) !=
" ")
goto L140;
5838 ctemp = crdbuf.substr(ibegin-1,lenbuf-ibegin);
5840 if (ierr > 0)
goto L180;
5851 scanf(crdbuf.c_str(),xk,stmp,uk,wk,a,b);
5854 if (k == 0)
goto L210;
5858 mnparm(k, cnamk, uk, wk, a, b, ierr);
5886 Double_urt a, f, s, t, y, s2, x2, x3, x4, y2, cz[3], xm, xy, x2y;
5896 for (i = 1; i <= 3; ++i) { cz[i-1] = 0; }
5898 if (npar2p < 3)
goto L10;
5902 for (i = 1; i <= npar2p; ++i) { xm += parx2p[i]; }
5911 for (i = 1; i <= npar2p; ++i) {
5923 a = (f*x4 - x2*x2)*x2 - f*(x3*x3);
5924 if (a == 0)
goto L10;
5925 cz[2] = (x2*(f*x2y - x2*y) - f*x3*xy) / a;
5926 cz[1] = (xy - x3*cz[2]) / x2;
5927 cz[0] = (y - x2*cz[2]) / f;
5928 if (npar2p == 3)
goto L6;
5929 sdev2p = y2 - (cz[0]*y + cz[1]*xy + cz[2]*x2y);
5930 if (sdev2p < 0) sdev2p = 0;
5933 cz[0] += xm*(xm*cz[2] - cz[1]);
5934 cz[1] -= xm*2*cz[2];
5936 for (i = 1; i <= 3; ++i) { coef2p[i] = cz[i-1]; }
5950 string chbuf2, chbufi;
5956 igo = m_userParameterFlag[i];
5961 yy = (pexti - alimi)*2 / (blimi - alimi) - 1;
5966 chbuf2 =
" IS AT ITS LOWER ALLOWED LIMIT.";
5969 chbuf2 =
" IS AT ITS UPPER ALLOWED LIMIT.";
5972 pexti = alimi + (blimi - alimi)*.5*(
URMath::Sin(a) + 1);
5974 if (yy2 > 1) chbuf2 =
" BROUGHT BACK INSIDE LIMITS.";
5975 ostringstream warning;
5976 warning <<
"VARIABLE " << i <<
" " << chbuf2;
5999 static string cdot =
".";
6000 static string cslash =
"/";
6003 Double_urt xmin, ymin, xmax, ymax, savx, savy, yprt;
6004 Double_urt bwidx, bwidy, xbest, ybest, ax, ay, bx, by;
6006 Int_urt iten, i, j, k, maxnx, maxny, iquit, ni, linodd;
6007 Int_urt nxbest, nybest, km1, ibk, isp1, nx, ny, ks, ix;
6008 string chmess, ctemp;
6016 if (maxnx < 10) maxnx = 10;
6018 if (maxny < 10) maxny = 10;
6019 if (nxypt <= 1)
return;
6025 for (i = 1; i <= km1; ++i) {
6028 for (j = 1; j <= ni; ++j) {
6029 if (ypt[j-1] > ypt[j])
continue;
6041 if (iquit == 0)
break;
6046 for (i = 1; i <= nxypt; ++i) {
6047 if (xpt[i-1] > xmax) xmax = xpt[i-1];
6048 if (xpt[i-1] < xmin) xmin = xpt[i-1];
6050 dxx = (xmax - xmin)*.001;
6053 mnbins(xmin, xmax, maxnx, xmin, xmax, nx, bwidx);
6055 ymin = ypt[nxypt-1];
6056 if (ymax == ymin) ymax = ymin + 1;
6057 dyy = (ymax - ymin)*.001;
6060 mnbins(ymin, ymax, maxny, ymin, ymax, ny, bwidy);
6063 if (chbest ==
' ')
goto L50;
6064 xbest = (xmax + xmin)*.5;
6065 ybest = (ymax + ymin)*.5;
6073 for (i = 1; i <= nxypt; ++i) {
6074 xpt[i-1] = ax*xpt[i-1] + bx;
6075 ypt[i-1] = any - ay*ypt[i-1] - by;
6077 nxbest =
Int_urt((ax*xbest + bx));
6078 nybest =
Int_urt((any - ay*ybest - by));
6085 for (i = 1; i <= ny; ++i) {
6086 for (ibk = 1; ibk <= nx; ++ibk) { cline[ibk-1] =
' '; }
6091 cline[nxbest-1] =
'.';
6092 if (i != 1 && i != nybest && i != ny)
goto L320;
6093 for (j = 1; j <= nx; ++j) { cline[j-1] =
'.'; }
6096 if (isp1 > nxypt)
goto L350;
6098 for (k = isp1; k <= nxypt; ++k) {
6100 if (ks > i)
goto L345;
6102 if (cline[ix-1] ==
'.')
goto L340;
6103 if (cline[ix-1] ==
' ')
goto L340;
6104 if (cline[ix-1] == chpt[k-1])
continue;
6111 cline[ix-1] = chpt[k-1];
6118 if (linodd == 1 || i == ny)
goto L380;
6121 std::printf(
" %s\n",ctemp.c_str());
6125 std::printf(
" %14.7g ..%s\n",yprt,ctemp.c_str());
6131 for (ibk = 1; ibk <= nx; ++ibk) {
6133 if (ibk % 10 == 1) cline[ibk-1] =
'/';
6135 std::printf(
" %s\n",cline);
6137 for (ibk = 1; ibk <= 12; ++ibk) {
6138 xvalus[ibk-1] = xmin +
Double_urt(ibk-1)*10*bwidx;
6141 iten = (nx + 9) / 10;
6142 for (ibk = 1; ibk <= iten; ++ibk) {
6143 std::printf(
" %9.4g", xvalus[ibk-1]);
6146 if (overpr) chmess =
" Overprint character is &";
6147 std::printf(
" ONE COLUMN=%13.7g%s\n",bwidx,chmess.c_str());
6177 if (iuext == 0)
goto L100;
6181 if (iint >
fNpar)
goto L100;
6187 if (iext == 0)
goto L100;
6188 if (iext >
fNu)
goto L100;
6190 iint = m_userParameterIdToInternalId[iext];
6195 nvl = m_userParameterFlag[iext];
6196 if (nvl < 0)
goto L100;
6199 chnam = m_userParameterName[iext];
6200 val = m_userParameterValue[iext];
6201 if (iint > 0) err =
fWerr[iint-1];
6203 xlolim =
fAlim[iext-1];
6204 xuplim =
fBlim[iext-1];
6210 chnam =
"undefined";
6233 static string cblank =
" ";
6234 static string cnambf =
" ";
6239 Int_urt nadd, i, k, l, m, ikode, ic, nc, ntrail, lbl;
6241 string colhdl[6], colhdu[6], cx2, cx3, cheval;
6245 *m_logStream <<
" THERE ARE CURRENTLY NO PARAMETERS DEFINED" << endl;
6251 ikode =
fISW[1] + 1;
6252 if (ikode > 3) ikode = 3;
6255 for (k = 1; k <= 6; ++k) {
6256 colhdu[k-1] =
"UNDEFINED";
6257 colhdl[k-1] =
"COLUMN HEAD";
6262 *m_logStream <<
" MINUIT TASK: " <<
fCtitl <<
'\n';
6265 if (fval ==
fUndefi) cheval =
" unknown ";
6284 *m_logStream <<
" FCN=" << cheval <<
" FROM " << setw(8) <<
fCfrom <<
" STATUS=" 6285 << setw(10) <<
fCstatu << setw(7) << nc <<
" " << setw(9) <<
fNfcn <<
" TOTAL\n" << endl;
6287 if (m == 0 || m == 2 ||
fDcovar == 0) {
6291 *m_logStream <<
" EDM=" << chedm <<
" STRATEGY=" << setw(2)
6298 *m_logStream <<
" EDM=" << chedm
6299 <<
" STRATEGY=" << setw(2) <<
fIstrat <<
" ERROR MATRIX UNCERTAINTY " 6300 << setw(5) << setprecision(1) << dc <<
" per cent\n";
6303 if (ikode == 0)
return;
6306 for (i = 1; i <=
fNu; ++i) {
6308 if (m_userParameterFlag[i] < 0)
continue;
6310 uint lastChar = m_userParameterName[i].size();
6311 uint icMax = (lastChar < 10 ) ? lastChar : 10;
6312 for (ic = icMax; ic >= 1; --ic) {
6314 if (m_userParameterName[i].substr(ic-1,1) !=
" ")
goto L16;
6319 if (lbl < ntrail) ntrail = lbl;
6321 nadd = ntrail / 2 + 1;
6324 colhdl[0] =
" ERROR ";
6325 colhdu[1] =
" PHYSICAL";
6326 colhdu[2] =
" LIMITS ";
6327 colhdl[1] =
" NEGATIVE ";
6328 colhdl[2] =
" POSITIVE ";
6332 colhdl[0] =
" ERROR ";
6333 colhdu[1] =
" INTERNAL ";
6334 colhdl[1] =
" STEP SIZE ";
6335 colhdu[2] =
" INTERNAL ";
6336 colhdl[2] =
" VALUE ";
6340 colhdl[0] =
" ERROR ";
6341 colhdu[1] =
" STEP ";
6342 colhdl[1] =
" SIZE ";
6343 colhdu[2] =
" FIRST ";
6344 colhdl[2] =
" DERIVATIVE ";
6347 colhdu[0] =
" PARABOLIC ";
6348 colhdl[0] =
" ERROR ";
6349 colhdu[1] =
" MINOS ";
6350 colhdu[2] =
"ERRORS ";
6351 colhdl[1] =
" NEGATIVE ";
6352 colhdl[2] =
" POSITIVE ";
6356 if (
fISW[1] < 3) colhdu[0] =
" APPROXIMATE ";
6357 if (
fISW[1] < 1) colhdu[0] =
" CURRENT GUESS";
6360 int maxNameSize = 0;
6361 for ( vector<string>::iterator name = m_userParameterName.begin();
6362 name != m_userParameterName.end();
6364 int nameSize = name->size();
6365 maxNameSize = max(maxNameSize,nameSize);
6367 maxNameSize = max(maxNameSize,9);
6370 *m_logStream <<
" EXT " << setw(maxNameSize) << left <<
"PARAMETER" <<
" " 6371 << setw(14) << colhdu[0].c_str() << setw(14) << colhdu[1].c_str() << setw(14) << colhdu[2].c_str()
6373 *m_logStream <<
" NO. " << setw(maxNameSize) << left <<
" NAME " <<
" VALUE " 6374 << setw(14) << colhdl[0].c_str() << setw(14) << colhdl[1].c_str() << setw(14) << colhdl[2].c_str()
6383 for (i = 1; i <=
fNu; ++i) {
6385 if (m_userParameterFlag[i] < 0)
continue;
6387 l = m_userParameterIdToInternalId[i];
6389 cnambf = cblank.substr(0,nadd) + m_userParameterName[i];
6390 if (l == 0)
goto L55;
6393 cx2 =
"PLEASE GET X..";
6394 cx3 =
"PLEASE GET X..";
6397 if (m_userParameterFlag[i] <= 1) {
6400 *m_logStream << setw(5) << right << i <<
" " << setw(maxNameSize) << m_userParameterName[i]
6401 << setw(14) << setprecision(5) << m_userParameterValue[i]
6402 << setw(14) << setprecision(5) << x1 << endl;
6418 cx3 =
"** at limit **";
6423 if (x2 == 0) cx2 =
" ";
6424 if (x2 ==
fUndefi) cx2 =
" at limit ";
6426 if (x3 == 0) cx3 =
" ";
6427 if (x3 ==
fUndefi) cx3 =
" at limit ";
6429 if (cx2 ==
"PLEASE GET X..") {
6431 ocx << setw(14) << setprecision(5) << x2;
6435 if (cx3 ==
"PLEASE GET X..") {
6437 ocx << setw(14) << setprecision(5) << x3;
6447 *m_logStream << setw(5) << right << i <<
" " << setw(maxNameSize) << m_userParameterName[i]
6448 << setw(14) << setprecision(5) << m_userParameterValue[i]
6449 << setw(14) << setprecision(5) << x1
6450 << cx2 << cx3 << endl;
6454 if (m_userParameterFlag[i] <= 1 || ikode == 3)
continue;
6457 *m_logStream <<
" WARNING - - ABOVE PARAMETER IS AT LIMIT.\n";
6463 colhdu[0] =
" constant ";
6465 if (m_userParameterFlag[i] > 0) colhdu[0] =
" fixed ";
6467 if (m_userParameterFlag[i] == 4 && ikode == 1) {
6474 *m_logStream << setw(5) << right << i <<
" " << setw(maxNameSize) << m_userParameterName[i]
6475 << setw(14) << setprecision(5) << m_userParameterValue[i]
6476 << setw(14) << colhdu[0].c_str()
6477 << setw(14) << setprecision(5) <<
fAlim[i-1]
6478 << setw(14) << setprecision(5) <<
fBlim[i-1]
6485 *m_logStream << setw(5) << right << i <<
" " << setw(maxNameSize) << m_userParameterName[i]
6486 << setw(14) << setprecision(5) << m_userParameterValue[i]
6487 << setw(14) << colhdu[0].c_str()
6493 std::printf(
" ERR DEF= %g\n",
fUp);
6494 *m_logStream <<
" ERR DEF= " <<
fUp;
6496 *m_logStream << endl;
6509 Double_urt dgmin, padd, pmin, pmax, dg, epspdf, epsmin;
6510 Int_urt ndex, i, j, ndexd, ip, ifault;
6511 string chbuff, ctemp;
6517 for (i = 1; i <=
fNpar; ++i) {
6518 ndex = i*(i + 1) / 2;
6519 if (
fVhmat[ndex-1] <= 0) {
6520 ostringstream warning;
6521 warning <<
"Negative diagonal element " << i <<
" in Error Matrix";
6527 dg = epspdf + 1 - dgmin;
6528 ostringstream warning;
6529 warning << dg <<
" added to diagonal of error matrix";
6535 for (i = 1; i <=
fNpar; ++i) {
6540 for (j = 1; j <= i; ++j) {
6549 for (ip = 2; ip <=
fNpar; ++ip) {
6554 if ( ( pmin <= 0 &&
fLwarn ) ||
fISW[4] >= 2) {
6555 std::printf(
" EIGENVALUES OF SECOND-DERIVATIVE MATRIX:\n");
6557 for (ip = 1; ip <=
fNpar; ++ip) {
6558 ostringstream localTemp;
6559 localTemp <<
fPstar[ip-1];
6560 ctemp += localTemp.str();
6562 std::printf(
"%s",ctemp.c_str());
6564 if (pmin > epspdf*pmax)
return;
6566 padd = pmax*.001 - pmin;
6567 for (ip = 1; ip <=
fNpar; ++ip) {
6568 ndex = ip*(ip + 1) / 2;
6569 fVhmat[ndex-1] *= padd + 1;
6572 ostringstream warning;
6573 warning <<
"MATRIX FORCED POS-DEF BY ADDING %f TO DIAGONAL. " << padd;
6595 for (i = 1; i <=
fNpar; ++i) {
fX[i-1] = pnew[i-1]; }
6603 for (j = 2; j <= nparp1; ++j) {
if (y[j-1] > y[jh-1]) jh = j; }
6604 fEDM = y[jh-1] - y[jl-1];
6605 if (
fEDM <= 0)
goto L45;
6606 for (i = 1; i <=
fNpar; ++i) {
6609 for (j = 2; j <= nparp1; ++j) {
6613 fDirin[i-1] = pbig - plit;
6618 std::printf(
" FUNCTION VALUE DOES NOT SEEM TO DEPEND ON ANY OF THE%d VARIABLE PARAMETERS.\n",
fNpar);
6619 std::printf(
" VERIFY THAT STEP SIZES ARE BIG ENOUGH AND CHECK FCN LOGIC.\n");
6620 std::printf(
" *******************************************************************************\n");
6621 std::printf(
" *******************************************************************************\n");
6642 if (val == 3)
goto L100;
6645 iseed = (iseed - k*53668)*40014 - k*12211;
6646 if (iseed < 0) iseed += 2147483563;
6678 for (i = 1; i <=
fNpar; ++i) {
6701 std::printf(
"mnsave is dummy in the base class URMinuit: Use URMinuitOld\n");
6716 Double_urt step, uhigh, xhreq, xlreq, ubest, fnext, unext, xh, xl;
6717 Int_urt ipar, iint, icall, ncall, nbins, nparx;
6718 Int_urt nxypt, nccall, iparwd;
6723 if (ncall <= 1) ncall = 41;
6724 if (ncall > 101) ncall = 101;
6730 iint = m_userParameterIdToInternalId[ipar];
6732 if (iparwd > 0)
goto L200;
6737 if (ipar >
fNu)
goto L900;
6739 iint = m_userParameterIdToInternalId[ipar];
6740 if (iint <= 0)
goto L100;
6744 ubest = m_userParameterValue[ipar];
6753 if (m_userParameterFlag[ipar] > 1)
goto L300;
6756 if (xlreq == xhreq)
goto L250;
6758 step = (xhreq - xlreq) /
Double_urt(ncall-1);
6761 xl = ubest -
fWerr[iint-1];
6762 xh = ubest + fWerr[iint-1];
6763 mnbins(xl, xh, ncall, unext, uhigh, nbins, step);
6768 if (xlreq == xhreq)
goto L350;
6773 if (xl >= xh)
goto L700;
6778 unext =
fAlim[ipar-1];
6782 for (icall = 1; icall <= nccall; ++icall) {
6784 m_userParameterValue[ipar] = unext;
6789 fXpt[nxypt-1] = unext;
6790 fYpt[nxypt-1] = fnext;
6791 fChpt[nxypt-1] =
'*';
6792 if (fnext <
fAmin) {
6803 m_userParameterValue[ipar] = ubest;
6807 std::printf(
"%dSCAN OF PARAMETER NO. %d, %s\n" 6808 ,
fNewpag,ipar,m_userParameterName[ipar].c_str());
6812 std::printf(
" REQUESTED RANGE OUTSIDE LIMITS FOR PARAMETER %d\n",ipar);
6814 if (iparwd <= 0)
goto L100;
6835 Double_urt dxdi, rnum, ftry, rnum1, rnum2, alpha;
6837 Int_urt ipar, iext, j, ifail, iseed, nparx, istep, ib, mxfail, mxstep;
6840 if (mxfail <= 0) mxfail =
fNpar*20 + 100;
6844 if (alpha <= 0) alpha = 3;
6846 std::printf(
" MNSEEK: MONTE CARLO MINIMIZATION USING METROPOLIS ALGORITHM\n");
6847 std::printf(
" TO STOP AFTER %6d SUCCESSIVE FAILURES, OR %7d STEPS\n",mxfail,mxstep);
6848 std::printf(
" MAXIMUM STEP SIZE IS %9.3f ERROR BARS.\n",alpha);
6860 for (ipar = 1; ipar <=
fNpar; ++ipar) {
6864 if (m_userParameterFlag[iext] > 1) {
6867 if (dxdi == 0) dxdi = 1;
6868 fDirin[ipar-1] = alpha*2*fWerr[ipar-1] / dxdi;
6870 fDirin[ipar-1] = 6.2831859999999997;
6877 for (istep = 1; istep <= mxstep; ++istep) {
6878 if (ifail >= mxfail)
break;
6879 for (ipar = 1; ipar <=
fNpar; ++ipar) {
6910 std::printf(
" MNSEEK: %5d SUCCESSIVE UNSUCCESSFUL TRIALS.\n",ifail);
6934 const char *cname[30] = {
6968 static string cprlev[5] = {
6969 "-1: NO OUTPUT EXCEPT FROM SHOW ",
6970 " 0: REDUCED OUTPUT ",
6971 " 1: NORMAL OUTPUT ",
6972 " 2: EXTRA OUTPUT FOR PROBLEM CASES",
6973 " 3: MAXIMUM OUTPUT "};
6975 static string cstrat[3] = {
6976 " 0: MINIMIZE THE NUMBER OF CALLS TO FUNCTION",
6977 " 1: TRY TO BALANCE SPEED AGAINST RELIABILITY",
6978 " 2: MAKE SURE MINIMUM TRUE, ERRORS CORRECT "};
6980 static string cdbopt[7] = {
6981 "REPORT ALL EXCEPTIONAL CONDITIONS ",
6982 "MNLINE: LINE SEARCH MINIMIZATION ",
6983 "MNDERI: FIRST DERIVATIVE CALCULATIONS ",
6984 "MNHESS: SECOND DERIVATIVE CALCULATIONS ",
6985 "MNMIGR: COVARIANCE MATRIX UPDATES ",
6986 "MNHES1: FIRST DERIVATIVE UNCERTAINTIES ",
6987 "MNCONT: MNCONTOUR PLOT (MNCROS SEARCH) "};
6994 Int_urt iset, iprm, i, jseed, kname, iseed, iunit, id, ii, kk;
6995 Int_urt ikseed, idbopt, igrain, iswsav, isw2;
6996 string cfname, cmode, ckind, cwarn, copt, ctemp, ctemp2;
6999 for (i = 1; i <= nntot; ++i) {
7001 ckind = ctemp.substr(0,3);
7002 ctemp2 =
fCword.substr(4,6);
7003 if (strstr(ctemp2.c_str(),ckind.c_str()))
goto L5;
7010 ctemp2 =
fCword.substr(0,3);
7011 if ( ctemp2.find(
"HEL") != string::npos)
goto L2000;
7012 if ( ctemp2.find(
"SHO") != string::npos)
goto L1000;
7013 if ( ctemp2.find(
"SET") == string::npos)
goto L1900;
7017 if (kname <= 0)
goto L1900;
7019 switch ((
int)kname) {
7037 case 18:
goto L3000;
7039 case 20:
goto L3000;
7044 case 25:
goto L3000;
7045 case 26:
goto L1900;
7055 if (iprm >
fNu)
goto L25;
7056 if (iprm <= 0)
goto L25;
7058 if (m_userParameterFlag[iprm] < 0)
goto L25;
7060 m_userParameterValue[iprm] =
fWord7[1];
7071 std::printf(
" UNDEFINED PARAMETER NUMBER. IGNORED.\n");
7102 for (i = 1; i <=
fNpar; ++i) {
7130 mnwarn(
"W",
"SHO",
"SHO");
7138 std::printf(
" MINUIT RANDOM NUMBER SEED SET TO %d\n",jseed);
7150 if (
fISW[4] > 0)
goto L1172;
7168 if (
fISW[4] >= 0)
goto L1220;
7173 if (
fISW[4] >= 0)
goto L1100;
7178 if (
fISW[4] >= 0)
goto L1100;
7189 if (idbopt > 6)
goto L288;
7191 fIdbg[idbopt] = iset;
7192 if (iset == 1)
fIdbg[0] = 1;
7195 for (
id = 0;
id <= 6; ++id) {
fIdbg[id] = iset; }
7198 mnwarn(
"D",
"SHO",
"SHO");
7201 std::printf(
" UNKNOWN DEBUG OPTION %d REQUESTED. IGNORED\n",idbopt);
7212 if (kname <= 0)
goto L1900;
7214 switch ((
int)kname) {
7224 case 10:
goto L1100;
7225 case 11:
goto L1110;
7226 case 12:
goto L1120;
7227 case 13:
goto L1130;
7228 case 14:
goto L1130;
7229 case 15:
goto L1150;
7230 case 16:
goto L1160;
7231 case 17:
goto L1170;
7232 case 18:
goto L1180;
7233 case 19:
goto L1190;
7234 case 20:
goto L1200;
7235 case 21:
goto L1210;
7236 case 22:
goto L1220;
7237 case 23:
goto L1100;
7238 case 24:
goto L1100;
7239 case 25:
goto L1250;
7240 case 26:
goto L1900;
7241 case 27:
goto L1270;
7242 case 28:
goto L1270;
7243 case 29:
goto L1290;
7244 case 30:
goto L1300;
7274 std::printf(
" ALLOWED PRINT LEVELS ARE:\n");
7275 std::printf(
" %s\n",cprlev[0].c_str());
7276 std::printf(
" %s\n",cprlev[1].c_str());
7277 std::printf(
" %s\n",cprlev[2].c_str());
7278 std::printf(
" %s\n",cprlev[3].c_str());
7279 std::printf(
" %s\n",cprlev[4].c_str());
7280 std::printf(
" CURRENT PRINTOUT LEVEL IS %s\n",cprlev[
fISW[4]+1].c_str());
7285 std::printf(
" NOGRAD IS SET. DERIVATIVES NOT COMPUTED IN FCN.\n");
7287 std::printf(
" GRAD IS SET. USER COMPUTES DERIVATIVES IN FCN.\n");
7292 std::printf(
" ERRORS CORRESPOND TO FUNCTION CHANGE OF %g\n",
fUp);
7316 cmode =
"BATCH MODE ";
7317 if (
fISW[5] == 1) cmode =
"INTERACTIVE MODE";
7318 if (! lname) cfname =
"unknown";
7319 std::printf(
" INPUT NOW BEING READ IN %s FROM UNIT NO. %d FILENAME: %s\n" 7320 ,cmode.c_str(),
fIsysrd,cfname.c_str());
7324 std::printf(
" PAGE WIDTH IS SET TO %d COLUMNS\n",
fNpagwd);
7328 std::printf(
" PAGE LENGTH IS SET TO %d LINES\n",
fNpagln);
7332 cwarn =
"SUPPRESSED";
7333 if (
fLwarn) cwarn =
"REPORTED ";
7334 std::printf(
"%s\n",cwarn.c_str());
7342 std::printf(
" MINUIT RNDM SEED IS CURRENTLY=%d\n",ikseed);
7349 std::printf(
" TITLE OF CURRENT TASK IS:%s\n",
fCtitl.c_str());
7353 std::printf(
" ALLOWED STRATEGIES ARE:\n");
7354 std::printf(
" %s\n",cstrat[0].c_str());
7355 std::printf(
" %s\n",cstrat[1].c_str());
7356 std::printf(
" %s\n",cstrat[2].c_str());
7358 std::printf(
" NOW USING STRATEGY %s\n",cstrat[
fIstrat].c_str());
7365 std::printf(
"%s\n",
fCovmes[0].c_str());
7373 std::printf(
" PAGE THROW CARRIAGE CONTROL = %d\n",
fNewpag);
7375 std::printf(
" NO PAGE THROWS IN MINUIT OUTPUT\n");
7380 for (ii = 1; ii <=
fNpar; ++ii) {
7381 if (
fErp[ii-1] > 0 ||
fErn[ii-1] < 0)
goto L1204;
7383 std::printf(
" THERE ARE NO MINOS ERRORS CURRENTLY VALID.\n");
7390 std::printf(
" FLOATING-POINT NUMBERS ASSUMED ACCURATE TO %g\n",
fEpsmac);
7394 std::printf(
" MINUIT PRIMARY OUTPUT TO UNIT %d\n",
fIsyswr);
7398 std::printf(
" THIS IS MINUIT VERSION:%s\n",
fCvrsn.c_str());
7402 for (
id = 0;
id <= 6; ++id) {
7404 if (
fIdbg[
id] >= 1) copt =
"ON ";
7405 std::printf(
" DEBUG OPTION %3d IS %3s :%s\n" 7406 ,
id,copt.c_str(),cdbopt[id].c_str());
7421 std::printf(
" THE COMMAND:%10s IS UNKNOWN.\n",
fCword.c_str());
7427 ctemp2 =
fCword.substr(3,7);
7428 if (strcmp(ctemp2.c_str(),
"SHO")) ckind =
"SHOW";
7430 std::printf(
" THE FORMAT OF THE %4s COMMAND IS:\n",ckind.c_str());
7431 std::printf(
" %s xxx [numerical arguments if any]\n",ckind.c_str());
7432 std::printf(
" WHERE xxx MAY BE ONE OF THE FOLLOWING:\n");
7433 for (kk = 1; kk <= nname; ++kk) {
7434 std::printf(
" %s\n",cname[kk-1]);
7440 std::printf(
" ABOVE COMMAND IS ILLEGAL. IGNORED\n");
7462 Double_urt dmin_, dxdi, yrho, f, ynpp1, aming, ypbar;
7463 Double_urt bestx, ystar, y1, y2, ystst, pb, wg;
7465 Int_urt npfn, i, j, k, jhold, ncycl, nparx;
7466 Int_urt nparp1, kg, jh, nf, jl, ns;
7468 if (
fNpar <= 0)
return;
7477 rho2 = rho1 + alpha*gamma;
7480 std::printf(
" START SIMPLEX MINIMIZATION. CONVERGENCE WHEN EDM .LT. %g\n",
fEpsi);
7482 for (i = 1; i <=
fNpar; ++i) {
7495 for (i = 1; i <=
fNpar; ++i) {
7507 if (f <= aming)
goto L6;
7509 if (kg == 1)
goto L8;
7513 if (nf < 3)
goto L4;
7523 if (ns < 6)
goto L4;
7527 if (aming < absmin) jl = i;
7528 if (aming < absmin) absmin = aming;
7548 for (i = 1; i <=
fNpar; ++i) {
7557 if (ystar >=
fAmin)
goto L70;
7559 for (i = 1; i <=
fNpar; ++i) {
7566 y1 = (ystar -
fSIMPy[jh-1])*rho2;
7567 y2 = (ystst -
fSIMPy[jh-1])*rho1;
7568 rho = (rho2*y1 - rho1*y2)*.5 / (y1 - y2);
7569 if (rho < rhomin)
goto L66;
7570 if (rho > rhomax) rho = rhomax;
7571 for (i = 1; i <=
fNpar; ++i) {
7577 if (yrho <
fSIMPy[jl-1] && yrho < ystst)
goto L65;
7578 if (ystst <
fSIMPy[jl-1])
goto L67;
7579 if (yrho >
fSIMPy[jl-1])
goto L66;
7585 if (ystst <
fSIMPy[jl-1])
goto L67;
7592 if (
fISW[4] < 2)
goto L50;
7593 if (
fISW[4] >= 3 || ncycl % 10 == 0) {
7599 if (ystar >=
fSIMPy[jh-1])
goto L73;
7602 if (jhold != jh)
goto L50;
7605 for (i = 1; i <=
fNpar; ++i) {
7611 if (ystst >
fSIMPy[jh-1])
goto L1;
7613 if (ystst <
fAmin)
goto L67;
7619 std::printf(
" SIMPLEX MINIMIZATION HAS CONVERGED.\n");
7625 std::printf(
" SIMPLEX TERMINATES WITHOUT CONVERGENCE.\n");
7631 for (i = 1; i <=
fNpar; ++i) {
7703 static string cpt =
" ABCDEFGHIJKLMNOPQRSTUVWXYZabcdefghijklmnopqrstuvwxyz1234567890./;:[]$%*_!@#&+()";
7706 l = strlen(cfname.c_str());
7707 for (i = 1; i <= l; ++i) {
7708 for (ic = 1; ic <= 80; ++ic) {
7709 if (cfname[i-1] == cpt[ic-1])
goto L100;
7741 if (n < 1)
goto L100;
7744 for (i = 1; i <= n; ++i) {
7746 if (si <= 0)
goto L100;
7749 for (i = 1; i <= n; ++i) {
7750 for (j = 1; j <= n; ++j) {
7755 for (i = 1; i <= n; ++i) {
7758 if (a[k + k*l] != 0)
fVERTq[k-1] = 1 / a[k + k*l];
7764 if (km1 < 0)
goto L100;
7765 else if (km1 == 0)
goto L50;
7768 for (j = 1; j <= km1; ++j) {
7774 if (k - n < 0)
goto L51;
7775 else if (k - n == 0)
goto L60;
7778 for (j = kp1; j <= n; ++j) {
7785 for (j = 1; j <= n; ++j) {
7786 for (k = j; k <= n; ++k) { a[j + k*l] +=
fVERTpp[j-1]*
fVERTq[k-1]; }
7790 for (j = 1; j <= n; ++j) {
7791 for (k = 1; k <= j; ++k) {
7793 a[j + k*l] = a[k + j*l];
7817 string copt = copt1;
7818 string corg = corg1;
7819 string cmes = cmes1;
7823 string englsh, ctyp;
7825 if (corg.substr(0,3) !=
"SHO" || cmes.substr(0,3) !=
"SHO") {
7831 std::printf(
" MINUIT WARNING IN %s\n",corg.c_str());
7832 std::printf(
" ============== %s\n",cmes.c_str());
7838 std::printf(
" MINUIT DEBUG FOR %s\n",corg.c_str());
7839 std::printf(
" =============== %s \n",cmes.c_str());
7864 englsh =
" WAS SUPPRESSED. ";
7865 if (
fNwrmes[ityp-1] > 1) englsh =
"S WERE SUPPRESSED.";
7866 std::printf(
" %5d MINUIT %s MESSAGE%s\n",
fNwrmes[ityp-1]
7867 ,ctyp.c_str(),englsh.c_str());
7871 std::printf(
" ONLY THE MOST RECENT 10 WILL BE LISTED BELOW.\n");
7875 std::printf(
" CALLS ORIGIN MESSAGE\n");
7876 for (i = 1; i <= nm; ++i) {
7878 if (ic > MAXMES) ic = 1;
7896 Int_urt ndex, ierr, i, j, k, l, ndiag, k1, iin;
7900 for (l = 1; l <=
fNpar; ++l) {
7901 ndex = l*(l + 1) / 2;
7905 if (m_userParameterFlag[i] > 1) {
7907 ba =
fBlim[i-1] - al;
7910 du1 = al + 0.5*(
URMath::Sin(
fX[l-1] + dx) + 1)*ba - m_userParameterValue[i];
7911 du2 = al + 0.5*(
URMath::Sin(
fX[l-1] - dx) + 1)*ba - m_userParameterValue[i];
7912 if (dx > 1) du1 = ba;
7920 for (i = 1; i <=
fNpar; ++i) {
7923 for (j = 1; j <= i; ++j) {
7931 for (iin = 1; iin <=
fNpar; ++iin) {
7932 ndiag = iin*(iin + 1) / 2;
7934 if (denom <= 1 && denom >= 0)
fGlobcc[iin-1] = 0;
static Double_urt ASin(Double_urt)
virtual Int_urt FixParameter(Int_urt parNo)
virtual void mnpars(const std::string &crdbuf, Int_urt &icondn)
static Double_urt Sin(Double_urt)
void mnwarn(const char *copt1, const char *corg1, const char *cmes1)
virtual void mnemat(Double_urt *emat, Int_urt ndim)
virtual void mnpout(Int_urt iuext, std::string &chnam, Double_urt &val, Double_urt &err, Double_urt &xlolim, Double_urt &xuplim, Int_urt &iuint) const
static Short_urt Max(Short_urt a, Short_urt b)
virtual Int_urt GetParameter(Int_urt parNo, Double_urt ¤tValue, Double_urt ¤tError) const
virtual Int_urt GetNumFreePars() const
virtual void mntiny(Double_urt epsp1, Double_urt &epsbak)
virtual void mnexcm(const std::string &command, Double_urt *plist, Int_urt llist, Int_urt &ierflg)
virtual void BuildArrays(Int_urt maxpar=15)
static Short_urt Min(Short_urt a, Short_urt b)
virtual Int_urt DefineParameter(Int_urt parNo, const std::string &name, Double_urt initVal, Double_urt initErr, Double_urt lowerLimit, Double_urt upperLimit)
Bool_urt mnunpt(const std::string &cfname)
virtual void mnexin(Double_urt *pint)
virtual void mnfixp(Int_urt iint, Int_urt &ierr)
virtual Int_urt SetPrintLevel(Int_urt printLevel=0)
static Double_urt Log10(Double_urt x)
virtual void mnvert(Double_urt *a, Int_urt l, Int_urt m, Int_urt n, Int_urt &ifail)
virtual void mnrn15(Double_urt &val, Int_urt &inseed)
virtual Int_urt GetNumFixedPars() const
static Double_urt ATan(Double_urt)
static Double_urt Log(Double_urt x)
virtual void mnpint(Double_urt &pexti, Int_urt i, Double_urt &pinti)
virtual void mnrset(Int_urt iopt)
virtual void mndxdi(Double_urt pint, Int_urt ipar, Double_urt &dxdi)
const Int_urt kDefaultMaximumInternalParameters
virtual Int_urt Release(Int_urt parNo)
virtual void mnrazz(Double_urt ynew, Double_urt *pnew, Double_urt *y, Int_urt &jh, Int_urt &jl)
virtual void mneig(Double_urt *a, Int_urt ndima, Int_urt n, Int_urt mits, Double_urt *work, Double_urt precis, Int_urt &ifault)
static Double_urt Sqrt(Double_urt x)
virtual void mnparm(Int_urt k, const std::string &cnamj, Double_urt uk, Double_urt wk, Double_urt a, Double_urt b, Int_urt &ierflg)
virtual void mncont(Int_urt ke1, Int_urt ke2, Int_urt nptu, Double_urt *xptu, Double_urt *yptu, Int_urt &ierrf)
virtual void mnfree(Int_urt k)
static Double_urt Power(Double_urt x, Double_urt y)
const string kUndefinedParameterName(")UNDEFINED")
virtual void mnline(Double_urt *start, Double_urt fstart, Double_urt *step, Double_urt slope, Double_urt toler)
std::string fWarmes[kMAXWARN]
virtual void mnerrs(Int_urt number, Double_urt &eplus, Double_urt &eminus, Double_urt &eparab, Double_urt &gcc)
static Short_urt Abs(Short_urt d)
virtual void mninit(Int_urt i1, Int_urt i2, Int_urt i3)
virtual Int_urt Command(const char *command)
virtual void mncros(Double_urt &aopt, Int_urt &iercr)
std::string fOrigin[kMAXWARN]
const Int_urt kDefaultMaximumExternalParameters
virtual void mninex(Double_urt *pint)
virtual Int_urt GetNumPars() const
virtual void mncntr(Int_urt ke1, Int_urt ke2, Int_urt &ierrf)
virtual void mnmatu(Int_urt kode)
virtual void DeleteArrays()
virtual void mnstat(Double_urt &fmin, Double_urt &fedm, Double_urt &errdef, Int_urt &npari, Int_urt &nparx, Int_urt &istat)
bool parameterFixed(Int_urt parameterNumber)
virtual void mncrck(std::string crdbuf, Int_urt maxcwd, std::string &comand, Int_urt &lnc, Int_urt mxp, Double_urt *plist, Int_urt &llist, Int_urt &ierr, Int_urt isyswr)
std::string fCfrom
Character to be plotted at the X,Y contour positions.
virtual void mnpfit(Double_urt *parx2p, Double_urt *pary2p, Int_urt npar2p, Double_urt *coef2p, Double_urt &sdev2p)
virtual void mnprin(Int_urt inkode, Double_urt fval)
virtual Int_urt SetErrorDef(Double_urt up)
virtual void mnmnot(Int_urt ilax, Int_urt ilax2, Double_urt &val2pl, Double_urt &val2mi)
virtual void SetMaxIterations(Int_urt maxiter=5000)
virtual Int_urt Eval(Int_urt npar, Double_urt *grad, Double_urt &fval, const std::vector< Double_urt > &par, Int_urt flag)
virtual void mnbins(Double_urt a1, Double_urt a2, Int_urt naa, Double_urt &bl, Double_urt &bh, Int_urt &nb, Double_urt &bwid)
virtual void mncalf(Double_urt *pvec, Double_urt &ycalf)
virtual void mneval(Double_urt anext, Double_urt &fnext, Int_urt &ierev)
virtual void mncomd(const std::string &crdbin, Int_urt &icondn)
virtual void SetFCN(URFcn *fcn)
virtual void mnplot(Double_urt *xpt, Double_urt *ypt, char *chpt, Int_urt nxypt, Int_urt npagwd, Int_urt npagln)
virtual void zeroPointers()
static Double_urt Cos(Double_urt)