31 const std::vector<double>& pdir,
double tlr,
unsigned int maxcalls)
const {
39 unsigned int npar = par.size();
40 unsigned int nfcn = 0;
43 double mgr_tlr = 0.5 * tlr;
52 unsigned int maxitr = 15;
54 double aminsv =
fFval;
55 double aim = aminsv + up;
59 std::vector<double> alsb(3, 0.), flsb(3, 0.);
65 std::cout<<
"MnFunctionCross for parameter "<<par.front()<<
"fmin = " << aminsv
66 <<
" contur value aim = (fmin + up) = " << aim << std::endl;
73 for(
unsigned int i = 0; i < par.size(); i++) {
74 unsigned int kex = par[i];
76 double zmid = pmid[i];
77 double zdir = pdir[i];
82 aulim =
std::min(aulim, (zlim-zmid)/zdir);
86 aulim =
std::min(aulim, (zlim-zmid)/zdir);
92 std::cout<<
"Largest allowed aulim "<< aulim << std::endl;
95 if(aulim < aopt+tla) limset =
true;
100 for(
unsigned int i = 0; i < npar; i++) {
102 std::cout <<
"MnFunctionCross: Set value for " << par[i] <<
" to " << pmid[i] << std::endl;
106 if (printLevel > 1) {
107 std::cout <<
"MnFunctionCross: parameter " << i <<
" set to " << pmid[i] << std::endl;
116 std::cout <<
"MnFunctionCross: after Migrad on n-1 minimum is " << min0 << std::endl;
122 if(limset ==
true && min0.
Fval() < aim)
127 flsb[0] = min0.
Fval();
128 flsb[0] =
std::max(flsb[0], aminsv + 0.1*up);
129 aopt =
sqrt(up/(flsb[0]-aminsv)) - 1.;
132 if(aopt > 1.) aopt = 1.;
133 if(aopt < -0.5) aopt = -0.5;
140 std::cout <<
"MnFunctionCross: flsb[0] = " << flsb[0] <<
" aopt = " << aopt << std::endl;
143 for(
unsigned int i = 0; i < npar; i++) {
145 std::cout <<
"MnFunctionCross: Set new value for " << par[i] <<
" from " << pmid[i] <<
" to " << pmid[i] + (aopt)*pdir[i] <<
" aopt = " << aopt << std::endl;
147 migrad.
SetValue(par[i], pmid[i] + (aopt)*pdir[i]);
149 if (printLevel > 1) {
150 std::cout <<
"MnFunctionCross: parameter " << i <<
" set to " << pmid[i] + (aopt)*pdir[i] << std::endl;
159 std::cout <<
"MnFunctionCross: after Migrad on n-1 minimum is " << min1 << std::endl;
165 if(limset ==
true && min1.
Fval() < aim)
170 flsb[1] = min1.
Fval();
171 double dfda = (flsb[1] - flsb[0])/(alsb[1] - alsb[0]);
174 std::cout <<
"aopt = " << aopt <<
" min1Val = " << flsb[1] <<
" dfda = " << dfda << std::endl;
182 std::cout <<
"MnFunctionCross: dfda < 0 - iterate from " << ipt <<
" to max of " << maxitr << std::endl;
186 unsigned int maxlk = maxitr - ipt;
187 for(
unsigned int it = 0; it < maxlk; it++) {
191 aopt = alsb[0] + 0.2*(it+1);
197 for(
unsigned int i = 0; i < npar; i++) {
199 std::cout <<
"MnFunctionCross: Set new value for " << par[i] <<
" to " << pmid[i] + (aopt)*pdir[i] <<
" aopt = " << aopt << std::endl;
201 migrad.
SetValue(par[i], pmid[i] + (aopt)*pdir[i]);
202 if (printLevel > 1) {
203 std::cout <<
"MnFunctionCross: parameter " << i <<
" set to " << pmid[i] + (aopt)*pdir[i] << std::endl;
207 min1 = migrad(maxcalls, mgr_tlr);
211 std::cout <<
"MnFunctionCross: after Migrad on n-1 minimum is " << min1 << std::endl;
212 std::cout <<
"nfcn = " << nfcn << std::endl;
218 if(limset ==
true && min1.
Fval() < aim)
222 flsb[1] = min1.
Fval();
223 dfda = (flsb[1] - flsb[0])/(alsb[1] - alsb[0]);
227 std::cout <<
"aopt = " << aopt <<
" min1Val = " << flsb[1] <<
" dfda = " << dfda << std::endl;
239 aopt = alsb[1] + (aim-flsb[1])/dfda;
242 std::cout <<
"MnFunctionCross: dfda > 0 : aopt = " << aopt << std::endl;
248 if(
fabs(aopt) > 1.) tla = tlr*
fabs(aopt);
249 if(adist < tla && fdist < tlf)
return MnCross(aopt, min1.
UserState(), nfcn);
251 double bmin =
std::min(alsb[0], alsb[1]) - 1.;
252 if(aopt < bmin) aopt = bmin;
253 double bmax =
std::max(alsb[0], alsb[1]) + 1.;
254 if(aopt > bmax) aopt = bmax;
262 for(
unsigned int i = 0; i < npar; i++) {
264 std::cout <<
"MnFunctionCross: Set new value for " << par[i] <<
" from " << pmid[i] <<
" to " << pmid[i] + (aopt)*pdir[i] <<
" aopt = " << aopt << std::endl;
266 migrad.
SetValue(par[i], pmid[i] + (aopt)*pdir[i]);
267 if (printLevel > 1) {
268 std::cout <<
"MnFunctionCross: parameter " << i <<
" set to " << pmid[i] + (aopt)*pdir[i] << std::endl;
276 std::cout <<
"MnFunctionCross: after Migrad on n-1 minimum is " << min2 << std::endl;
277 std::cout <<
"nfcn = " << nfcn << std::endl;
283 if(limset ==
true && min2.
Fval() < aim)
288 flsb[2] = min2.
Fval();
292 double ecarmn =
fabs(flsb[2] - aim);
294 unsigned int ibest = 2;
295 unsigned int iworst = 0;
296 unsigned int noless = 0;
298 for(
unsigned int i = 0; i < 3; i++) {
299 double ecart =
fabs(flsb[i] - aim);
308 if(flsb[i] < aim) noless++;
312 std::cout <<
"MnFunctionCross: have three points : nless < aim = " << noless <<
" ibest = " << ibest <<
" iworst = " << iworst << std::endl;
318 if(noless == 1 || noless == 2)
goto L500;
323 if(noless == 3 && ibest != 2) {
327 std::cout <<
"MnFunctionCross: all three points below - look again fir positive slope " << std::endl;
334 flsb[iworst] = flsb[2];
335 alsb[iworst] = alsb[2];
336 dfda = (flsb[1] - flsb[0])/(alsb[1] - alsb[0]);
338 std::cout <<
"MnFunctionCross: new straight line using point 1-2 - dfda = " << dfda << std::endl;
352 std::cout <<
"MnFunctionCross: parabola fit: iteration " << ipt << std::endl;
355 double coeff1 = parbol.
C();
356 double coeff2 = parbol.
B();
357 double coeff3 = parbol.
A();
358 double determ = coeff2*coeff2 - 4.*coeff3*(coeff1 - aim);
361 std::cout <<
"MnFunctionCross: parabola fit: a = " << coeff3 <<
" b = "
362 << coeff2 <<
" c = " << coeff1 <<
" determ = " << determ << std::endl;
366 double rt =
sqrt(determ);
367 double x1 = (-coeff2 + rt)/(2.*coeff3);
368 double x2 = (-coeff2 - rt)/(2.*coeff3);
369 double s1 = coeff2 + 2.*x1*coeff3;
370 double s2 = coeff2 + 2.*x2*coeff3;
373 std::cout <<
"MnFunctionCross: parabola fit: x1 = " << x1 <<
" x2 = "
374 << x2 <<
" s1 = " << s1 <<
" s2 = " << s2 << std::endl;
378 if(s1*s2 > 0.)
MN_INFO_MSG(
"MnFunctionCross problem 1");
388 std::cout <<
"MnFunctionCross: parabola fit: aopt = " << aopt <<
" slope = "
389 << slope << std::endl;
394 if(
fabs(aopt) > 1.) tla = tlr*
fabs(aopt);
397 std::cout <<
"MnFunctionCross: Delta(aopt) = " <<
fabs(aopt - alsb[ibest]) <<
" tla = "
398 << tla <<
"Delta(F) = " <<
fabs(flsb[ibest] - aim) <<
" tlf = " << tlf << std::endl;
402 if(
fabs(aopt - alsb[ibest]) < tla &&
fabs(flsb[ibest] - aim) < tlf)
410 unsigned int ileft = 3;
411 unsigned int iright = 3;
412 unsigned int iout = 3;
415 ecarmn =
fabs(aim-flsb[0]);
416 for(
unsigned int i = 0; i < 3; i++) {
417 double ecart =
fabs(flsb[i] - aim);
422 if(ecart > ecarmx) ecarmx = ecart;
424 if(iright == 3) iright = i;
425 else if(flsb[i] > flsb[iright]) iout = i;
430 }
else if(ileft == 3) ileft = i;
431 else if(flsb[i] < flsb[ileft]) iout = i;
439 std::cout <<
"MnFunctionCross: ileft = " << ileft <<
" iright = "
440 << iright <<
" iout = " << iout <<
" ibest = " << ibest << std::endl;
445 if(ecarmx > 10.*
fabs(flsb[iout] - aim))
446 aopt = 0.5*(aopt + 0.5*(alsb[iright] + alsb[ileft]));
449 double smalla = 0.1*tla;
450 if(slope*smalla > tlf) smalla = tlf/slope;
451 double aleft = alsb[ileft] + smalla;
452 double aright = alsb[iright] - smalla;
455 if(aopt < aleft) aopt = aleft;
456 if(aopt > aright) aopt = aright;
457 if(aleft > aright) aopt = 0.5*(aleft + aright);
467 for(
unsigned int i = 0; i < npar; i++) {
469 std::cout <<
"MnFunctionCross: Set new value for " << par[i] <<
" from " << pmid[i] <<
" to " << pmid[i] + (aopt)*pdir[i] <<
" aopt = " << aopt << std::endl;
471 migrad.
SetValue(par[i], pmid[i] + (aopt)*pdir[i]);
472 if (printLevel > 1) {
473 std::cout <<
"MnFunctionCross: parameter " << i <<
" set to " << pmid[i] + (aopt)*pdir[i] << std::endl;
476 min2 = migrad(maxcalls, mgr_tlr);
479 std::cout <<
"MnFunctionCross: after Migrad on n-1 minimum is " << min2 << std::endl;
480 std::cout <<
"nfcn = " << nfcn << std::endl;
486 if(limset ==
true && min2.
Fval() < aim)
492 flsb[iout] = min2.
Fval();
494 }
while(ipt < maxitr);
virtual double Up() const =0
Error definition of the function.
API class for minimization using Variable Metric technology ("MIGRAD"); allows for user interaction: ...
static Vc_ALWAYS_INLINE int_v min(const int_v &x, const int_v &y)
double A() const
Accessor to the coefficient of the quadratic term.
const MnStrategy & fStrategy
const MnUserParameterState & UserState() const
determines the relative floating point arithmetic precision.
double LowerLimit() const
static const double x2[5]
double C() const
Accessor to the coefficient of the constant term.
class holding the full result of the minimization; both internal and external (MnUserParameterState) ...
This class defines a parabola of the form a*x*x + b*x + c.
const MnUserParameterState & fState
double UpperLimit() const
MnCross operator()(const std::vector< unsigned int > &, const std::vector< double > &, const std::vector< double > &, double, unsigned int) const
bool HasUpperLimit() const
VecExpr< UnaryOp< Fabs< T >, VecExpr< A, T, D >, T >, T, D > fabs(const VecExpr< A, T, D > &rhs)
double B() const
Accessor to the coefficient of the linear term.
const MnMachinePrecision & Precision() const
unsigned int Strategy() const
const MinuitParameter & Parameter(unsigned int i) const
static const double x1[5]
double Eps() const
eps returns the smallest possible number so that 1.+eps > 1.
bool HasReachedCallLimit() const
static Vc_ALWAYS_INLINE int_v max(const int_v &x, const int_v &y)
void SetValue(unsigned int, double)
bool HasLowerLimit() const
API class for defining three levels of strategies: low (0), medium (1), high (>=2); acts on: Migrad (...