35 #if !defined(R__SOLARIS) && !defined(R__ACC) && !defined(R__FBSD)
67 if(
x==0.0)
return 0.0;
69 return log(
x+ax*
sqrt(1.+1./(ax*ax)));
80 if(
x==0.0)
return 0.0;
82 return log(
x+ax*
sqrt(1.-1./(ax*ax)));
93 return log((1+
x)/(1-
x))/2;
118 const Double_t c[20] = {0.42996693560813697, 0.40975987533077106,
119 -0.01858843665014592, 0.00145751084062268,-0.00014304184442340,
120 0.00001588415541880,-0.00000190784959387, 0.00000024195180854,
121 -0.00000003193341274, 0.00000000434545063,-0.00000000060578480,
122 0.00000000008612098,-0.00000000001244332, 0.00000000000182256,
123 -0.00000000000027007, 0.00000000000004042,-0.00000000000000610,
124 0.00000000000000093,-0.00000000000000014, 0.00000000000000002};
127 t=
h=
y=
s=
a=alfa=b1=b2=b0=0.;
131 }
else if (
x == -1) {
140 a = -pi3+hf*(b1*b1-b2*b2);
146 }
else if (t <= -0.5) {
170 for (
Int_t i=19;i>=0;i--){
171 b0 =
c[i] + alfa*b1-b2;
175 h = -(
s*(b0-
h*b2)+
a);
207 Double_t kConst = 0.8862269254527579;
212 Double_t erfi, derfi, y0,y1,dy0,dy1;
217 for (
Int_t iter=0; iter<kMaxit; iter++) {
220 if (
TMath::Abs(dy1) < kEps) {
if (
x < 0)
return -erfi;
else return erfi;}
225 if(
TMath::Abs(derfi/erfi) < kEps) {
if (
x < 0)
return -erfi;
else return erfi;}
249 if (
n <= 0)
return 1.;
268 const Double_t w2 = 1.41421356237309505;
270 const Double_t p10 = 2.4266795523053175e+2, q10 = 2.1505887586986120e+2,
271 p11 = 2.1979261618294152e+1, q11 = 9.1164905404514901e+1,
272 p12 = 6.9963834886191355e+0, q12 = 1.5082797630407787e+1,
273 p13 =-3.5609843701815385e-2, q13 = 1;
275 const Double_t p20 = 3.00459261020161601e+2, q20 = 3.00459260956983293e+2,
276 p21 = 4.51918953711872942e+2, q21 = 7.90950925327898027e+2,
277 p22 = 3.39320816734343687e+2, q22 = 9.31354094850609621e+2,
278 p23 = 1.52989285046940404e+2, q23 = 6.38980264465631167e+2,
279 p24 = 4.31622272220567353e+1, q24 = 2.77585444743987643e+2,
280 p25 = 7.21175825088309366e+0, q25 = 7.70001529352294730e+1,
281 p26 = 5.64195517478973971e-1, q26 = 1.27827273196294235e+1,
282 p27 =-1.36864857382716707e-7, q27 = 1;
284 const Double_t p30 =-2.99610707703542174e-3, q30 = 1.06209230528467918e-2,
285 p31 =-4.94730910623250734e-2, q31 = 1.91308926107829841e-1,
286 p32 =-2.26956593539686930e-1, q32 = 1.05167510706793207e+0,
287 p33 =-2.78661308609647788e-1, q33 = 1.98733201817135256e+0,
288 p34 =-2.23192459734184686e-2, q34 = 1;
339 if (
x > 0)
return 0.5 +0.5*
h;
381 if (
a <= 0 ||
x <= 0)
return 0;
389 for (
Int_t i=1; i<=itmax; i++) {
393 if (
Abs(
d) < fpmin)
d = fpmin;
395 if (
Abs(
c) < fpmin)
c = fpmin;
399 if (
Abs(del-1) < eps)
break;
417 if (
a <= 0 ||
x <= 0)
return 0;
450 if (
sigma == 0)
return 1.e30;
453 if (arg < -39.0 || arg > 39.0)
return 0.0;
455 if (!norm)
return res;
456 return res/(2.50662827463100024*
sigma);
471 if (
sigma <= 0)
return 0;
473 if (!norm)
return den;
520 if( av0 >= av1 && av0 >= av2 ) {
526 else if (av1 >= av0 && av1 >= av2) {
542 Double_t d = amax *
Sqrt(1.+foofrac*foofrac+barfrac*barfrac);
616 if (ndf <= 0)
return 0;
619 if (chi2 < 0)
return 0;
669 }
else if (u < 0.755) {
672 }
else if (u < 6.8116) {
678 for (
Int_t j=0; j<maxj;j++) {
681 p = 2*(
r[0] -
r[1] +
r[2] -
r[3]);
791 if (!
a || !
b || na <= 2 || nb <= 2) {
792 Error(
"KolmogorovTest",
"Sets must have more than 2 points");
808 for (
Int_t i=0;i<na+nb;i++) {
812 if (ia >= na) {ok =
kTRUE;
break;}
813 }
else if (
a[ia] >
b[ib]) {
816 if (ib >= nb) {ok =
kTRUE;
break;}
820 while(ia < na &&
a[ia] ==
x) {
824 while(ib < nb &&
b[ib] ==
x) {
828 if (ia >= na) {ok =
kTRUE;
break;}
829 if (ib >= nb) {ok =
kTRUE;
break;}
843 printf(
" Kolmogorov Probability = %g, Max Dist = %g\n",prob,rdmax);
877 if ((
sigma < 0 || lg < 0) || (
sigma==0 && lg==0)) {
882 return lg * 0.159154943 / (xx*xx + lg*lg /4);
890 x = xx /
sigma / 1.41421356;
891 y = lg / 2 /
sigma / 1.41421356;
910 Double_t c[6] = { 1.0117281, -0.75197147, 0.012557727, 0.010022008, -0.00024206814, 0.00000050084806};
911 Double_t s[6] = { 1.393237, 0.23115241, -0.15535147, 0.0062183662, 0.000091908299, -0.00000062752596};
912 Double_t t[6] = { 0.31424038, 0.94778839, 1.5976826, 2.2795071, 3.0206370, 3.8897249};
919 Double_t xlim0, xlim1, xlim2, xlim3, xlim4;
920 Double_t a0=0, d0=0, d2=0, e0=0, e2=0, e4=0, h0=0, h2=0, h4=0, h6=0;
921 Double_t p0=0, p2=0, p4=0, p6=0, p8=0, z0=0, z2=0, z4=0, z6=0, z8=0;
922 Double_t xp[6], xm[6], yp[6], ym[6];
923 Double_t mq[6], pq[6], mf[6], pf[6];
938 xlim3 = 3.097 *
y - 0.45;
940 xlim4 = 18.1 *
y + 1.65;
949 k = yrrtpi / (xq + yq);
950 }
else if ( abx > xlim1 ) {
957 d = rrtpi / (d0 + xq*(d2 + xq));
958 k =
d *
y * (a0 + xq);
959 }
else if ( abx > xlim2 ) {
962 h0 = 0.5625 + yq * (4.5 + yq * (10.5 + yq * (6.0 + yq)));
964 h2 = -4.5 + yq * (9.0 + yq * ( 6.0 + yq * 4.0));
965 h4 = 10.5 - yq * (6.0 - yq * 6.0);
966 h6 = -6.0 + yq * 4.0;
967 e0 = 1.875 + yq * (8.25 + yq * (5.5 + yq));
968 e2 = 5.25 + yq * (1.0 + yq * 3.0);
971 d = rrtpi / (h0 + xq * (h2 + xq * (h4 + xq * (h6 + xq))));
972 k =
d *
y * (e0 + xq * (e2 + xq * (e4 + xq)));
973 }
else if ( abx < xlim3 ) {
976 z0 = 272.1014 +
y * (1280.829 +
y *
986 z2 = 211.678 +
y * (902.3066 +
y *
994 z4 = 78.86585 +
y * (308.1852 +
y *
998 (80.39278 +
y * 10.0)
1000 z6 = 22.03523 +
y * (55.02933 +
y *
1002 (53.59518 +
y * 10.0)
1004 z8 = 1.496460 +
y * (13.39880 +
y * 5.0);
1005 p0 = 153.5168 +
y * (549.3954 +
y *
1012 (4.264678 +
y * 0.3183291)
1014 p2 = -34.16955 +
y * (-1.322256+
y *
1019 (12.79458 +
y * 1.2733163)
1021 p4 = 2.584042 +
y * (10.46332 +
y *
1024 (12.79568 +
y * 1.9099744)
1026 p6 = -0.07272979 +
y * (0.9377051 +
y *
1027 (4.266322 +
y * 1.273316));
1028 p8 = 0.0005480304 +
y * 0.3183291;
1030 d = 1.7724538 / (z0 + xq * (z2 + xq * (z4 + xq * (z6 + xq * (z8 + xq)))));
1031 k =
d * (p0 + xq * (p2 + xq * (p4 + xq * (p6 + xq * p8))));
1034 ypy0q = ypy0 * ypy0;
1036 for (j = 0; j <= 5; j++) {
1039 mf[j] = 1.0 / (mq[j] + ypy0q);
1041 ym[j] = mf[j] * ypy0;
1044 pf[j] = 1.0 / (pq[j] + ypy0q);
1046 yp[j] = pf[j] * ypy0;
1048 if ( abx <= xlim4 ) {
1049 for (j = 0; j <= 5; j++) {
1050 k = k +
c[j]*(ym[j]+yp[j]) -
s[j]*(xm[j]-xp[j]) ;
1054 for ( j = 0; j <= 5; j++) {
1056 (mq[j] * mf[j] - y0 * ym[j])
1057 +
s[j] * yf * xm[j]) / (mq[j]+y0q)
1058 + (
c[j] * (pq[j] * pf[j] - y0 * yp[j])
1059 -
s[j] * yf * xp[j]) / (pq[j]+y0q);
1061 k =
y * k +
exp( -xq );
1064 return k / 2.506628 /
sigma;
1087 Double_t r,
s,t,p,
q,
d,ps3,ps33,qs2,u,
v,tmp,lnu,lnv,su,sv,y1,y2,y3;
1091 if (coef[3] == 0)
return complex;
1092 r = coef[2]/coef[3];
1093 s = coef[1]/coef[3];
1094 t = coef[0]/coef[3];
1097 q = (2*
r*
r*
r)/27.0 - (
r*
s)/3 + t;
1186 if (type<1 || type>9){
1187 printf(
"illegal value of type\n");
1193 if (index) ind = index;
1196 isAllocated =
kTRUE;
1205 for (
Int_t i=0; i<nprob; i++){
1215 nppm =
n*prob[i]-0.5;
1239 if (
type == 4) {
a = 0;
b = 1; }
1240 else if (
type == 5) {
a = 0.5;
b = 0.5; }
1241 else if (
type == 6) {
a = 0.;
b = 0.; }
1242 else if (
type == 7) {
a = 1.;
b = 1.; }
1243 else if (
type == 8) {
a = 1./3.;
b =
a; }
1244 else if (
type == 9) {
a = 3./8.;
b =
a; }
1248 nppm =
a + prob[i] * (
n + 1 -
a -
b);
1256 int first = (j > 0 && j <=
n) ? j-1 : ( j <=0 ) ? 0 :
n-1;
1257 int second = (j > 0 && j <
n) ? j : ( j <=0 ) ? 0 :
n-1;
1291 if (Narr <= 0)
return;
1292 double *localArr1 =
new double[Narr];
1293 int *localArr2 =
new int[Narr];
1297 for(iEl = 0; iEl < Narr; iEl++) {
1298 localArr1[iEl] = arr1[iEl];
1299 localArr2[iEl] = iEl;
1302 for (iEl = 0; iEl < Narr; iEl++) {
1303 for (iEl2 = Narr-1; iEl2 > iEl; --iEl2) {
1304 if (localArr1[iEl2-1] < localArr1[iEl2]) {
1305 double tmp = localArr1[iEl2-1];
1306 localArr1[iEl2-1] = localArr1[iEl2];
1307 localArr1[iEl2] = tmp;
1309 int tmp2 = localArr2[iEl2-1];
1310 localArr2[iEl2-1] = localArr2[iEl2];
1311 localArr2[iEl2] = tmp2;
1316 for (iEl = 0; iEl < Narr; iEl++) {
1317 arr2[iEl] = localArr2[iEl];
1319 delete [] localArr2;
1320 delete [] localArr1;
1330 if (Narr <= 0)
return;
1331 double *localArr1 =
new double[Narr];
1332 int *localArr2 =
new int[Narr];
1336 for (iEl = 0; iEl < Narr; iEl++) {
1337 localArr1[iEl] = arr1[iEl];
1338 localArr2[iEl] = iEl;
1341 for (iEl = 0; iEl < Narr; iEl++) {
1342 for (iEl2 = Narr-1; iEl2 > iEl; --iEl2) {
1343 if (localArr1[iEl2-1] > localArr1[iEl2]) {
1344 double tmp = localArr1[iEl2-1];
1345 localArr1[iEl2-1] = localArr1[iEl2];
1346 localArr1[iEl2] = tmp;
1348 int tmp2 = localArr2[iEl2-1];
1349 localArr2[iEl2-1] = localArr2[iEl2];
1350 localArr2[iEl2] = tmp2;
1355 for (iEl = 0; iEl < Narr; iEl++) {
1356 arr2[iEl] = localArr2[iEl];
1358 delete [] localArr2;
1359 delete [] localArr1;
1404 const Double_t p1=1.0, p2=3.5156229, p3=3.0899424,
1405 p4=1.2067492, p5=0.2659732, p6=3.60768e-2, p7=4.5813e-3;
1407 const Double_t q1= 0.39894228, q2= 1.328592e-2, q3= 2.25319e-3,
1408 q4=-1.57565e-3, q5= 9.16281e-3, q6=-2.057706e-2,
1409 q7= 2.635537e-2, q8=-1.647633e-2, q9= 3.92377e-3;
1419 result = p1+
y*(p2+
y*(p3+
y*(p4+
y*(p5+
y*(p6+
y*p7)))));
1438 const Double_t p1=-0.57721566, p2=0.42278420, p3=0.23069756,
1439 p4= 3.488590e-2, p5=2.62698e-3, p6=1.0750e-4, p7=7.4e-6;
1441 const Double_t q1= 1.25331414, q2=-7.832358e-2, q3= 2.189568e-2,
1442 q4=-1.062446e-2, q5= 5.87872e-3, q6=-2.51540e-3, q7=5.3208e-4;
1445 Error(
"TMath::BesselK0",
"*K0* Invalid argument x = %g\n",
x);
1456 result = (
exp(-
x)/
sqrt(
x))*(q1+
y*(q2+
y*(q3+
y*(q4+
y*(q5+
y*(q6+
y*q7))))));
1472 const Double_t p1=0.5, p2=0.87890594, p3=0.51498869,
1473 p4=0.15084934, p5=2.658733e-2, p6=3.01532e-3, p7=3.2411e-4;
1475 const Double_t q1= 0.39894228, q2=-3.988024e-2, q3=-3.62018e-3,
1476 q4= 1.63801e-3, q5=-1.031555e-2, q6= 2.282967e-2,
1477 q7=-2.895312e-2, q8= 1.787654e-2, q9=-4.20059e-3;
1487 result =
x*(p1+
y*(p2+
y*(p3+
y*(p4+
y*(p5+
y*(p6+
y*p7))))));
1490 result = (
exp(ax)/
sqrt(ax))*(q1+
y*(q2+
y*(q3+
y*(q4+
y*(q5+
y*(q6+
y*(q7+
y*(q8+
y*q9))))))));
1491 if (
x < 0) result = -result;
1507 const Double_t p1= 1., p2= 0.15443144, p3=-0.67278579,
1508 p4=-0.18156897, p5=-1.919402e-2, p6=-1.10404e-3, p7=-4.686e-5;
1510 const Double_t q1= 1.25331414, q2= 0.23498619, q3=-3.655620e-2,
1511 q4= 1.504268e-2, q5=-7.80353e-3, q6= 3.25614e-3, q7=-6.8245e-4;
1514 Error(
"TMath::BesselK1",
"*K1* Invalid argument x = %g\n",
x);
1525 result = (
exp(-
x)/
sqrt(
x))*(q1+
y*(q2+
y*(q3+
y*(q4+
y*(q5+
y*(q6+
y*q7))))));
1538 if (
x <= 0 ||
n < 0) {
1539 Error(
"TMath::BesselK",
"*K* Invalid argument(s) (n,x) = (%d, %g)\n",
n,
x);
1551 for (
Int_t j=1; j<
n; j++) {
1568 const Double_t kBigPositive = 1.e10;
1569 const Double_t kBigNegative = 1.e-10;
1572 Error(
"TMath::BesselI",
"*I* Invalid argument (n,x) = (%d, %g)\n",
n,
x);
1579 if (
x == 0)
return 0;
1587 for (
Int_t j=
m; j>=1; j--) {
1593 result *= kBigNegative;
1595 bip *= kBigNegative;
1597 if (j==
n) result=bip;
1601 if ((
x < 0) && (
n%2 == 1)) result = -result;
1613 const Double_t p1 = 57568490574.0, p2 = -13362590354.0, p3 = 651619640.7;
1614 const Double_t p4 = -11214424.18, p5 = 77392.33017, p6 = -184.9052456;
1615 const Double_t p7 = 57568490411.0, p8 = 1029532985.0, p9 = 9494680.718;
1616 const Double_t p10 = 59272.64853, p11 = 267.8532712;
1619 const Double_t q2 = -0.1098628627e-2, q3 = 0.2734510407e-4;
1620 const Double_t q4 = -0.2073370639e-5, q5 = 0.2093887211e-6;
1621 const Double_t q6 = -0.1562499995e-1, q7 = 0.1430488765e-3;
1622 const Double_t q8 = -0.6911147651e-5, q9 = 0.7621095161e-6;
1623 const Double_t q10 = 0.934935152e-7, q11 = 0.636619772;
1625 if ((ax=
fabs(
x)) < 8) {
1627 result1 = p1 +
y*(p2 +
y*(p3 +
y*(p4 +
y*(p5 +
y*p6))));
1628 result2 = p7 +
y*(p8 +
y*(p9 +
y*(p10 +
y*(p11 +
y))));
1629 result = result1/result2;
1634 result1 = 1 +
y*(q2 +
y*(q3 +
y*(q4 +
y*q5)));
1635 result2 = q6 +
y*(q7 +
y*(q8 +
y*(q9 -
y*q10)));
1636 result =
sqrt(q11/ax)*(
cos(xx)*result1-z*
sin(xx)*result2);
1648 const Double_t p1 = 72362614232.0, p2 = -7895059235.0, p3 = 242396853.1;
1649 const Double_t p4 = -2972611.439, p5 = 15704.48260, p6 = -30.16036606;
1650 const Double_t p7 = 144725228442.0, p8 = 2300535178.0, p9 = 18583304.74;
1651 const Double_t p10 = 99447.43394, p11 = 376.9991397;
1654 const Double_t q2 = 0.183105e-2, q3 = -0.3516396496e-4;
1655 const Double_t q4 = 0.2457520174e-5, q5 = -0.240337019e-6;
1656 const Double_t q6 = 0.04687499995, q7 = -0.2002690873e-3;
1657 const Double_t q8 = 0.8449199096e-5, q9 = -0.88228987e-6;
1658 const Double_t q10 = 0.105787412e-6, q11 = 0.636619772;
1660 if ((ax=
fabs(
x)) < 8) {
1662 result1 =
x*(p1 +
y*(p2 +
y*(p3 +
y*(p4 +
y*(p5 +
y*p6)))));
1663 result2 = p7 +
y*(p8 +
y*(p9 +
y*(p10 +
y*(p11 +
y))));
1664 result = result1/result2;
1669 result1 = 1 +
y*(q2 +
y*(q3 +
y*(q4 +
y*q5)));
1670 result2 = q6 +
y*(q7 +
y*(q8 +
y*(q9 +
y*q10)));
1671 result =
sqrt(q11/ax)*(
cos(xx)*result1-z*
sin(xx)*result2);
1672 if (
x < 0) result = -result;
1683 const Double_t p1 = -2957821389., p2 = 7062834065.0, p3 = -512359803.6;
1684 const Double_t p4 = 10879881.29, p5 = -86327.92757, p6 = 228.4622733;
1685 const Double_t p7 = 40076544269., p8 = 745249964.8, p9 = 7189466.438;
1686 const Double_t p10 = 47447.26470, p11 = 226.1030244, p12 = 0.636619772;
1689 const Double_t q2 = -0.1098628627e-2, q3 = 0.2734510407e-4;
1690 const Double_t q4 = -0.2073370639e-5, q5 = 0.2093887211e-6;
1691 const Double_t q6 = -0.1562499995e-1, q7 = 0.1430488765e-3;
1692 const Double_t q8 = -0.6911147651e-5, q9 = 0.7621095161e-6;
1693 const Double_t q10 = -0.934945152e-7, q11 = 0.636619772;
1697 result1 = p1 +
y*(p2 +
y*(p3 +
y*(p4 +
y*(p5 +
y*p6))));
1698 result2 = p7 +
y*(p8 +
y*(p9 +
y*(p10 +
y*(p11 +
y))));
1704 result1 = 1 +
y*(q2 +
y*(q3 +
y*(q4 +
y*q5)));
1705 result2 = q6 +
y*(q7 +
y*(q8 +
y*(q9 +
y*q10)));
1706 result =
sqrt(q11/
x)*(
sin(xx)*result1+z*
cos(xx)*result2);
1717 const Double_t p1 = -0.4900604943e13, p2 = 0.1275274390e13;
1718 const Double_t p3 = -0.5153438139e11, p4 = 0.7349264551e9;
1719 const Double_t p5 = -0.4237922726e7, p6 = 0.8511937935e4;
1720 const Double_t p7 = 0.2499580570e14, p8 = 0.4244419664e12;
1721 const Double_t p9 = 0.3733650367e10, p10 = 0.2245904002e8;
1722 const Double_t p11 = 0.1020426050e6, p12 = 0.3549632885e3;
1725 const Double_t q2 = 0.183105e-2, q3 = -0.3516396496e-4;
1726 const Double_t q4 = 0.2457520174e-5, q5 = -0.240337019e-6;
1727 const Double_t q6 = 0.04687499995, q7 = -0.2002690873e-3;
1728 const Double_t q8 = 0.8449199096e-5, q9 = -0.88228987e-6;
1729 const Double_t q10 = 0.105787412e-6, q11 = 0.636619772;
1733 result1 =
x*(p1 +
y*(p2 +
y*(p3 +
y*(p4 +
y*(p5 +
y*p6)))));
1734 result2 = p7 +
y*(p8 +
y*(p9 +
y*(p10 +
y*(p11 +
y*(p12+
y)))));
1740 result1 = 1 +
y*(q2 +
y*(q3 +
y*(q4 +
y*q5)));
1741 result2 = q6 +
y*(q7 +
y*(q8 +
y*(q9 +
y*q10)));
1742 result =
sqrt(q11/
x)*(
sin(xx)*result1+z*
cos(xx)*result2);
1754 const Int_t n1 = 15;
1755 const Int_t n2 = 25;
1756 const Double_t c1[16] = { 1.00215845609911981, -1.63969292681309147,
1757 1.50236939618292819, -.72485115302121872,
1758 .18955327371093136, -.03067052022988,
1759 .00337561447375194, -2.6965014312602e-4,
1760 1.637461692612e-5, -7.8244408508e-7,
1761 3.021593188e-8, -9.6326645e-10,
1762 2.579337e-11, -5.8854e-13,
1763 1.158e-14, -2
e-16 };
1764 const Double_t c2[26] = { .99283727576423943, -.00696891281138625,
1765 1.8205103787037e-4, -1.063258252844e-5,
1766 9.8198294287e-7, -1.2250645445e-7,
1767 1.894083312e-8, -3.44358226e-9,
1768 7.1119102e-10, -1.6288744e-10,
1769 4.065681e-11, -1.091505e-11,
1770 3.12005e-12, -9.4202e-13,
1771 2.9848e-13, -9.872e-14,
1772 3.394e-14, -1.208e-14,
1773 4.44e-15, -1.68e-15,
1792 for (i = n1; i >= 0; --i) {
1793 b0 =
c1[i] + alfa*b1 - b2;
1805 for (i = n2; i >= 0; --i) {
1806 b0 =
c2[i] + alfa*b1 - b2;
1823 const Int_t n3 = 16;
1824 const Int_t n4 = 22;
1825 const Double_t c3[17] = { .5578891446481605, -.11188325726569816,
1826 -.16337958125200939, .32256932072405902,
1827 -.14581632367244242, .03292677399374035,
1828 -.00460372142093573, 4.434706163314e-4,
1829 -3.142099529341e-5, 1.7123719938e-6,
1830 -7.416987005e-8, 2.61837671e-9,
1831 -7.685839e-11, 1.9067e-12,
1832 -4.052e-14, 7.5e-16,
1834 const Double_t c4[23] = { 1.00757647293865641, .00750316051248257,
1835 -7.043933264519e-5, 2.66205393382e-6,
1836 -1.8841157753e-7, 1.949014958e-8,
1837 -2.6126199e-9, 4.236269e-10,
1838 -7.955156e-11, 1.679973e-11,
1839 -3.9072e-12, 9.8543e-13,
1840 -2.6636e-13, 7.645e-14,
1841 -2.313e-14, 7.33e-15,
1844 -4
e-17, 2
e-17,-1
e-17 };
1855 }
else if (
v <= 0.3) {
1860 for (i = 1; i <= i1; ++i) {
1861 h = -
h*
y / ((2*i+ 1)*(2*i + 3));
1871 for (i = n3; i >= 0; --i) {
1872 b0 =
c3[i] + alfa*b1 - b2;
1883 for (i = n4; i >= 0; --i) {
1884 b0 = c4[i] + alfa*b1 - b2;
1910 for (i=1; i<=60;i++) {
1911 r*=(
x/(2*i+1))*(
x/(2*i+1));
1919 for (i=1; i<=
km; i++) {
1920 r*=(2*i-1)*(2*i-1)/
x/
x;
1927 for (i=1; i<=16; i++) {
1928 r=0.125*
r*(2.0*i-1.0)*(2.0*i-1.0)/(i*
x);
1934 sl0=-2.0/(
pi*
x)*
s+bi0;
1952 for (i=1; i<=60;i++) {
1953 r*=
x*
x/(4.0*i*i-1.0);
1962 for (i=1; i<=
km; i++) {
1963 r*=(2*i+3)*(2*i+1)/
x/
x;
1967 sl1=2.0/
pi*(-1.0+1.0/(
x*
x)+3.0*
s/(
x*
x*
x*
x));
1971 for (i=1; i<=16; i++) {
1972 r=-0.125*
r*(4.0-(2.0*i-1.0)*(2.0*i-1.0))/(i*
x);
2006 d = 1.0 - qab*
x/qap;
2010 for (
m=1;
m<=itmax;
m++) {
2019 aa = -(
a+
m)*(qab +
m)*
x/((
a+
m2)*(qap+
m2));
2030 Info(
"TMath::BetaCf",
"a or b too big, or itmax too small, a=%g, b=%g, x=%g, h=%g, itmax=%d",
2046 if ((
x<0) || (
x>1) || (p<=0) || (
q<=0)){
2047 Error(
"TMath::BetaDist",
"parameter value outside allowed range");
2064 if ((
x<0) || (
x>1) || (p<=0) || (
q<=0)){
2065 Error(
"TMath::BetaDistI",
"parameter value outside allowed range");
2086 if (k==0 ||
n==k)
return 1;
2110 if(k <= 0)
return 1.0;
2111 if(k >
n)
return 0.0;
2159 Double_t c[]={0, 0.01, 0.222222, 0.32, 0.4, 1.24, 2.2,
2160 4.67, 6.66, 6.73, 13.32, 60.0, 70.0,
2161 84.0, 105.0, 120.0, 127.0, 140.0, 175.0,
2162 210.0, 252.0, 264.0, 294.0, 346.0, 420.0,
2163 462.0, 606.0, 672.0, 707.0, 735.0, 889.0,
2164 932.0, 966.0, 1141.0, 1182.0, 1278.0, 1740.0,
2172 if (ndf <= 0)
return 0;
2185 if (ch >
c[6]*ndf + 6)
2192 p1 = 1 + ch * (
c[7]+ch);
2193 p2 = ch * (
c[9] + ch * (
c[8] + ch));
2194 t = -0.5 + (
c[7] + 2 * ch) / p1 - (
c[9] + ch * (
c[10] + 3 * ch)) / p2;
2195 ch = ch - (1 -
TMath::Exp(
a +
g + 0.5 * ch + cp * aa) *p2 / p1) / t;
2200 if (ch <
e)
return ch;
2203 for (
Int_t i=0; i<maxit; i++){
2210 a = 0.5 * t -
b * cp;
2211 s1 = (
c[19] +
a * (
c[17] +
a * (
c[14] +
a * (
c[13] +
a * (
c[12] +
c[11] *
a))))) /
c[24];
2212 s2 = (
c[24] +
a * (
c[29] +
a * (
c[32] +
a * (
c[33] +
c[35] *
a)))) /
c[37];
2213 s3 = (
c[19] +
a * (
c[25] +
a * (
c[28] +
c[31] *
a))) /
c[37];
2214 s4 = (
c[20] +
a * (
c[27] +
c[34] *
a) + cp * (
c[22] +
a * (
c[30] +
c[36] *
a))) /
c[38];
2215 s5 = (
c[13] +
c[21] *
a + cp * (
c[18] +
c[26] *
a)) /
c[37];
2216 s6 = (
c[15] + cp * (
c[23] +
c[16] * cp)) /
c[38];
2217 ch = ch + t * (1 + 0.5 * t *
s1 -
b * cp * (
s1 -
b * (s2 -
b * (s3 -
b * (s4 -
b * (s5 -
b * s6))))));
2311 Error(
"TMath::GammaDist",
"illegal parameter values x = %f , gamma = %f beta = %f",
x,
gamma,
beta);
2399 if ((
x<theta) || (
sigma<=0) || (
m<=0)) {
2400 Error(
"TMath::Lognormal",
"illegal parameter values");
2418 if ((p<=0)||(p>=1)) {
2419 Error(
"TMath::NormQuantile",
"probability outside (0, 1)");
2423 Double_t a0 = 3.3871328727963666080e0;
2424 Double_t a1 = 1.3314166789178437745e+2;
2425 Double_t a2 = 1.9715909503065514427e+3;
2426 Double_t a3 = 1.3731693765509461125e+4;
2427 Double_t a4 = 4.5921953931549871457e+4;
2428 Double_t a5 = 6.7265770927008700853e+4;
2429 Double_t a6 = 3.3430575583588128105e+4;
2430 Double_t a7 = 2.5090809287301226727e+3;
2431 Double_t b1 = 4.2313330701600911252e+1;
2432 Double_t b2 = 6.8718700749205790830e+2;
2433 Double_t b3 = 5.3941960214247511077e+3;
2434 Double_t b4 = 2.1213794301586595867e+4;
2435 Double_t b5 = 3.9307895800092710610e+4;
2436 Double_t b6 = 2.8729085735721942674e+4;
2437 Double_t b7 = 5.2264952788528545610e+3;
2438 Double_t c0 = 1.42343711074968357734e0;
2442 Double_t c4 = 1.27045825245236838258e0;
2443 Double_t c5 = 2.41780725177450611770e-1;
2444 Double_t c6 = 2.27238449892691845833e-2;
2445 Double_t c7 = 7.74545014278341407640e-4;
2446 Double_t d1 = 2.05319162663775882187e0;
2447 Double_t d2 = 1.67638483018380384940e0;
2448 Double_t d3 = 6.89767334985100004550e-1;
2449 Double_t d4 = 1.48103976427480074590e-1;
2450 Double_t d5 = 1.51986665636164571966e-2;
2451 Double_t d6 = 5.47593808499534494600e-4;
2452 Double_t d7 = 1.05075007164441684324e-9;
2453 Double_t e0 = 6.65790464350110377720e0;
2454 Double_t e1 = 5.46378491116411436990e0;
2455 Double_t e2 = 1.78482653991729133580e0;
2456 Double_t e3 = 2.96560571828504891230e-1;
2457 Double_t e4 = 2.65321895265761230930e-2;
2458 Double_t e5 = 1.24266094738807843860e-3;
2459 Double_t e6 = 2.71155556874348757815e-5;
2460 Double_t e7 = 2.01033439929228813265e-7;
2462 Double_t f2 = 1.36929880922735805310e-1;
2463 Double_t f3 = 1.48753612908506148525e-2;
2464 Double_t f4 = 7.86869131145613259100e-4;
2465 Double_t f5 = 1.84631831751005468180e-5;
2466 Double_t f6 = 1.42151175831644588870e-7;
2467 Double_t f7 = 2.04426310338993978564e-15;
2478 quantile =
q* (((((((a7 *
r + a6) *
r + a5) *
r + a4) *
r + a3)
2479 *
r + a2) *
r + a1) *
r + a0) /
2480 (((((((b7 *
r + b6) *
r + b5) *
r + b4) *
r + b3)
2481 *
r + b2) *
r + b1) *
r + 1.);
2492 quantile=(((((((c7 *
r + c6) *
r + c5) *
r + c4) *
r +
c3)
2493 *
r +
c2) *
r +
c1) *
r + c0) /
2494 (((((((d7 *
r + d6) *
r + d5) *
r + d4) *
r + d3)
2495 *
r + d2) *
r + d1) *
r + 1);
2498 quantile=(((((((e7 *
r + e6) *
r + e5) *
r + e4) *
r + e3)
2499 *
r + e2) *
r + e1) *
r + e0) /
2500 (((((((f7 *
r + f6) *
r + f5) *
r + f4) *
r + f3)
2501 *
r + f2) *
r +
f1) *
r + 1);
2503 if (
q<0) quantile=-quantile;
2523 for(i=
n-2; i>-1; i--) {
2530 if(i1==-1)
return kFALSE;
2534 for(i=
n-1;i>i1;i--) {
2545 for(i=0;i<(
n-i1-1)/2;i++) {
2639 if (ndf<1 || p>=1 || p<=0) {
2640 Error(
"TMath::StudentQuantile",
"illegal parameter values");
2643 if ((lower_tail && p>0.5)||(!lower_tail && p<0.5)){
2645 q=2*(lower_tail ? (1-p) : p);
2648 q=2*(lower_tail? p : (1-p));
2668 if (ndf<5)
c+=0.3*(ndf-4.5)*(
x+0.6);
2669 c+=(((0.05*
d*
x-5.)*
x-7.)*
x-2.)*
x +
b;
2670 y=(((((0.4*
y+6.3)*
y+36.)*
y+94.5)/
c -
y-3.)/
b+1)*
x;
2675 y=((1./(((ndf+6.)/(ndf*
y)-0.089*
d-0.822)*(ndf+2.)*3)+0.5/(ndf+4.))*
y-1.)*
2676 (ndf+1.)/(ndf+2.)+1/
y;
2681 if(neg) quantile=-quantile;
2777 if (
x < ac[0])
v = 0;
2778 else if (
x >=ac[8])
v = 1;
2781 k =
Int_t(xx*ac[10]);
2782 v =
TMath::Min(wcm[k] + (xx - k*ac[9])*(wcm[k+1]-wcm[k])*ac[10], 1.);
2808 Double_t BKMNX1 = 0.02, BKMNY1 = 0.05, BKMNX2 = 0.12, BKMNY2 = 0.05,
2809 BKMNX3 = 0.22, BKMNY3 = 0.05, BKMXX1 = 0.1 , BKMXY1 = 1,
2810 BKMXX2 = 0.2 , BKMXY2 = 1 , BKMXX3 = 0.3 , BKMXY3 = 1;
2812 Double_t FBKX1 = 2/(BKMXX1-BKMNX1), FBKX2 = 2/(BKMXX2-BKMNX2),
2813 FBKX3 = 2/(BKMXX3-BKMNX3), FBKY1 = 2/(BKMXY1-BKMNY1),
2814 FBKY2 = 2/(BKMXY2-BKMNY2), FBKY3 = 2/(BKMXY3-BKMNY3);
2816 Double_t FNINV[] = {0, 1, 0.5, 0.33333333, 0.25, 0.2};
2818 Double_t EDGEC[]= {0, 0, 0.16666667e+0, 0.41666667e-1, 0.83333333e-2,
2819 0.13888889e-1, 0.69444444e-2, 0.77160493e-3};
2821 Double_t U1[] = {0, 0.25850868e+0, 0.32477982e-1, -0.59020496e-2,
2822 0. , 0.24880692e-1, 0.47404356e-2,
2823 -0.74445130e-3, 0.73225731e-2, 0. ,
2824 0.11668284e-2, 0. , -0.15727318e-2,-0.11210142e-2};
2826 Double_t U2[] = {0, 0.43142611e+0, 0.40797543e-1, -0.91490215e-2,
2827 0. , 0.42127077e-1, 0.73167928e-2,
2828 -0.14026047e-2, 0.16195241e-1, 0.24714789e-2,
2829 0.20751278e-2, 0. , -0.25141668e-2,-0.14064022e-2};
2831 Double_t U3[] = {0, 0.25225955e+0, 0.64820468e-1, -0.23615759e-1,
2832 0. , 0.23834176e-1, 0.21624675e-2,
2833 -0.26865597e-2, -0.54891384e-2, 0.39800522e-2,
2834 0.48447456e-2, -0.89439554e-2, -0.62756944e-2,-0.24655436e-2};
2836 Double_t U4[] = {0, 0.12593231e+1, -0.20374501e+0, 0.95055662e-1,
2837 -0.20771531e-1, -0.46865180e-1, -0.77222986e-2,
2838 0.32241039e-2, 0.89882920e-2, -0.67167236e-2,
2839 -0.13049241e-1, 0.18786468e-1, 0.14484097e-1};
2841 Double_t U5[] = {0, -0.24864376e-1, -0.10368495e-2, 0.14330117e-2,
2842 0.20052730e-3, 0.18751903e-2, 0.12668869e-2,
2843 0.48736023e-3, 0.34850854e-2, 0. ,
2844 -0.36597173e-3, 0.19372124e-2, 0.70761825e-3, 0.46898375e-3};
2846 Double_t U6[] = {0, 0.35855696e-1, -0.27542114e-1, 0.12631023e-1,
2847 -0.30188807e-2, -0.84479939e-3, 0. ,
2848 0.45675843e-3, -0.69836141e-2, 0.39876546e-2,
2849 -0.36055679e-2, 0. , 0.15298434e-2, 0.19247256e-2};
2851 Double_t U7[] = {0, 0.10234691e+2, -0.35619655e+1, 0.69387764e+0,
2852 -0.14047599e+0, -0.19952390e+1, -0.45679694e+0,
2853 0. , 0.50505298e+0};
2854 Double_t U8[] = {0, 0.21487518e+2, -0.11825253e+2, 0.43133087e+1,
2855 -0.14500543e+1, -0.34343169e+1, -0.11063164e+1,
2856 -0.21000819e+0, 0.17891643e+1, -0.89601916e+0,
2857 0.39120793e+0, 0.73410606e+0, 0. ,-0.32454506e+0};
2859 Double_t V1[] = {0, 0.27827257e+0, -0.14227603e-2, 0.24848327e-2,
2860 0. , 0.45091424e-1, 0.80559636e-2,
2861 -0.38974523e-2, 0. , -0.30634124e-2,
2862 0.75633702e-3, 0.54730726e-2, 0.19792507e-2};
2864 Double_t V2[] = {0, 0.41421789e+0, -0.30061649e-1, 0.52249697e-2,
2865 0. , 0.12693873e+0, 0.22999801e-1,
2866 -0.86792801e-2, 0.31875584e-1, -0.61757928e-2,
2867 0. , 0.19716857e-1, 0.32596742e-2};
2869 Double_t V3[] = {0, 0.20191056e+0, -0.46831422e-1, 0.96777473e-2,
2870 -0.17995317e-2, 0.53921588e-1, 0.35068740e-2,
2871 -0.12621494e-1, -0.54996531e-2, -0.90029985e-2,
2872 0.34958743e-2, 0.18513506e-1, 0.68332334e-2,-0.12940502e-2};
2874 Double_t V4[] = {0, 0.13206081e+1, 0.10036618e+0, -0.22015201e-1,
2875 0.61667091e-2, -0.14986093e+0, -0.12720568e-1,
2876 0.24972042e-1, -0.97751962e-2, 0.26087455e-1,
2877 -0.11399062e-1, -0.48282515e-1, -0.98552378e-2};
2879 Double_t V5[] = {0, 0.16435243e-1, 0.36051400e-1, 0.23036520e-2,
2880 -0.61666343e-3, -0.10775802e-1, 0.51476061e-2,
2881 0.56856517e-2, -0.13438433e-1, 0. ,
2882 0. , -0.25421507e-2, 0.20169108e-2,-0.15144931e-2};
2884 Double_t V6[] = {0, 0.33432405e-1, 0.60583916e-2, -0.23381379e-2,
2885 0.83846081e-3, -0.13346861e-1, -0.17402116e-2,
2886 0.21052496e-2, 0.15528195e-2, 0.21900670e-2,
2887 -0.13202847e-2, -0.45124157e-2, -0.15629454e-2, 0.22499176e-3};
2889 Double_t V7[] = {0, 0.54529572e+1, -0.90906096e+0, 0.86122438e-1,
2890 0. , -0.12218009e+1, -0.32324120e+0,
2891 -0.27373591e-1, 0.12173464e+0, 0. ,
2892 0. , 0.40917471e-1};
2894 Double_t V8[] = {0, 0.93841352e+1, -0.16276904e+1, 0.16571423e+0,
2895 0. , -0.18160479e+1, -0.50919193e+0,
2896 -0.51384654e-1, 0.21413992e+0, 0. ,
2897 0. , 0.66596366e-1};
2899 Double_t W1[] = {0, 0.29712951e+0, 0.97572934e-2, 0. ,
2900 -0.15291686e-2, 0.35707399e-1, 0.96221631e-2,
2901 -0.18402821e-2, -0.49821585e-2, 0.18831112e-2,
2902 0.43541673e-2, 0.20301312e-2, -0.18723311e-2,-0.73403108e-3};
2904 Double_t W2[] = {0, 0.40882635e+0, 0.14474912e-1, 0.25023704e-2,
2905 -0.37707379e-2, 0.18719727e+0, 0.56954987e-1,
2906 0. , 0.23020158e-1, 0.50574313e-2,
2907 0.94550140e-2, 0.19300232e-1};
2909 Double_t W3[] = {0, 0.16861629e+0, 0. , 0.36317285e-2,
2910 -0.43657818e-2, 0.30144338e-1, 0.13891826e-1,
2911 -0.58030495e-2, -0.38717547e-2, 0.85359607e-2,
2912 0.14507659e-1, 0.82387775e-2, -0.10116105e-1,-0.55135670e-2};
2914 Double_t W4[] = {0, 0.13493891e+1, -0.26863185e-2, -0.35216040e-2,
2915 0.24434909e-1, -0.83447911e-1, -0.48061360e-1,
2916 0.76473951e-2, 0.24494430e-1, -0.16209200e-1,
2917 -0.37768479e-1, -0.47890063e-1, 0.17778596e-1, 0.13179324e-1};
2919 Double_t W5[] = {0, 0.10264945e+0, 0.32738857e-1, 0. ,
2920 0.43608779e-2, -0.43097757e-1, -0.22647176e-2,
2921 0.94531290e-2, -0.12442571e-1, -0.32283517e-2,
2922 -0.75640352e-2, -0.88293329e-2, 0.52537299e-2, 0.13340546e-2};
2924 Double_t W6[] = {0, 0.29568177e-1, -0.16300060e-2, -0.21119745e-3,
2925 0.23599053e-2, -0.48515387e-2, -0.40797531e-2,
2926 0.40403265e-3, 0.18200105e-2, -0.14346306e-2,
2927 -0.39165276e-2, -0.37432073e-2, 0.19950380e-2, 0.12222675e-2};
2929 Double_t W8[] = {0, 0.66184645e+1, -0.73866379e+0, 0.44693973e-1,
2930 0. , -0.14540925e+1, -0.39529833e+0,
2931 -0.44293243e-1, 0.88741049e-1};
2934 if (rkappa <0.01 || rkappa >12) {
2935 Error(
"Vavilov distribution",
"illegal value of kappa");
2943 Double_t x,
y, xx, yy,
x2,
x3, y2, y3,
xy, p2, p3, q2, q3, pq;
2944 if (rkappa >= 0.29) {
2949 AC[0] = (-0.032227*beta2-0.074275)*rkappa + (0.24533*beta2+0.070152)*wk + (-0.55610*beta2-3.1579);
2950 AC[8] = (-0.013483*beta2-0.048801)*rkappa + (-1.6921*beta2+8.3656)*wk + (-0.73275*beta2-3.5226);
2953 for (j=1; j<=4; j++) {
2954 DRK[j+1] = DRK[1]*DRK[j];
2955 DSIGM[j+1] = DSIGM[1]*DSIGM[j];
2956 ALFA[j+1] = (FNINV[j]-beta2*FNINV[j+1])*DRK[j];
2960 HC[2]=ALFA[3]*DSIGM[3];
2961 HC[3]=(3*ALFA[2]*ALFA[2] + ALFA[4])*DSIGM[4]-3;
2962 HC[4]=(10*ALFA[2]*ALFA[3]+ALFA[5])*DSIGM[5]-10*
HC[2];
2966 for (j=2; j<=7; j++)
2968 HC[8]=0.39894228*
HC[1];
2970 else if (rkappa >=0.22) {
2973 x = 1+(rkappa-BKMXX3)*FBKX3;
2987 AC[1] = W1[1] + W1[2]*
x + W1[4]*
x3 + W1[5]*
y + W1[6]*y2 + W1[7]*y3 +
2988 W1[8]*
xy + W1[9]*p2 + W1[10]*p3 + W1[11]*q2 + W1[12]*q3 + W1[13]*pq;
2989 AC[2] = W2[1] + W2[2]*
x + W2[3]*
x2 + W2[4]*
x3 + W2[5]*
y + W2[6]*y2 +
2990 W2[8]*
xy + W2[9]*p2 + W2[10]*p3 + W2[11]*q2;
2991 AC[3] = W3[1] + W3[3]*
x2 + W3[4]*
x3 + W3[5]*
y + W3[6]*y2 + W3[7]*y3 +
2992 W3[8]*
xy + W3[9]*p2 + W3[10]*p3 + W3[11]*q2 + W3[12]*q3 + W3[13]*pq;
2993 AC[4] = W4[1] + W4[2]*
x + W4[3]*
x2 + W4[4]*
x3 + W4[5]*
y + W4[6]*y2 + W4[7]*y3 +
2994 W4[8]*
xy + W4[9]*p2 + W4[10]*p3 + W4[11]*q2 + W4[12]*q3 + W4[13]*pq;
2995 AC[5] = W5[1] + W5[2]*
x + W5[4]*
x3 + W5[5]*
y + W5[6]*y2 + W5[7]*y3 +
2996 W5[8]*
xy + W5[9]*p2 + W5[10]*p3 + W5[11]*q2 + W5[12]*q3 + W5[13]*pq;
2997 AC[6] = W6[1] + W6[2]*
x + W6[3]*
x2 + W6[4]*
x3 + W6[5]*
y + W6[6]*y2 + W6[7]*y3 +
2998 W6[8]*
xy + W6[9]*p2 + W6[10]*p3 + W6[11]*q2 + W6[12]*q3 + W6[13]*pq;
2999 AC[8] = W8[1] + W8[2]*
x + W8[3]*
x2 + W8[5]*
y + W8[6]*y2 + W8[7]*y3 + W8[8]*
xy;
3001 }
else if (rkappa >= 0.12) {
3004 x = 1 + (rkappa-BKMXX2)*FBKX2;
3018 AC[1] = V1[1] + V1[2]*
x + V1[3]*
x2 + V1[5]*
y + V1[6]*y2 + V1[7]*y3 +
3019 V1[9]*p2 + V1[10]*p3 + V1[11]*q2 + V1[12]*q3;
3020 AC[2] = V2[1] + V2[2]*
x + V2[3]*
x2 + V2[5]*
y + V2[6]*y2 + V2[7]*y3 +
3021 V2[8]*
xy + V2[9]*p2 + V2[11]*q2 + V2[12]*q3;
3022 AC[3] = V3[1] + V3[2]*
x + V3[3]*
x2 + V3[4]*
x3 + V3[5]*
y + V3[6]*y2 + V3[7]*y3 +
3023 V3[8]*
xy + V3[9]*p2 + V3[10]*p3 + V3[11]*q2 + V3[12]*q3 + V3[13]*pq;
3024 AC[4] = V4[1] + V4[2]*
x + V4[3]*
x2 + V4[4]*
x3 + V4[5]*
y + V4[6]*y2 + V4[7]*y3 +
3025 V4[8]*
xy + V4[9]*p2 + V4[10]*p3 + V4[11]*q2 + V4[12]*q3;
3026 AC[5] = V5[1] + V5[2]*
x + V5[3]*
x2 + V5[4]*
x3 + V5[5]*
y + V5[6]*y2 + V5[7]*y3 +
3027 V5[8]*
xy + V5[11]*q2 + V5[12]*q3 + V5[13]*pq;
3028 AC[6] = V6[1] + V6[2]*
x + V6[3]*
x2 + V6[4]*
x3 + V6[5]*
y + V6[6]*y2 + V6[7]*y3 +
3029 V6[8]*
xy + V6[9]*p2 + V6[10]*p3 + V6[11]*q2 + V6[12]*q3 + V6[13]*pq;
3030 AC[7] = V7[1] + V7[2]*
x + V7[3]*
x2 + V7[5]*
y + V7[6]*y2 + V7[7]*y3 +
3031 V7[8]*
xy + V7[11]*q2;
3032 AC[8] = V8[1] + V8[2]*
x + V8[3]*
x2 + V8[5]*
y + V8[6]*y2 + V8[7]*y3 +
3033 V8[8]*
xy + V8[11]*q2;
3037 if (rkappa >=0.02) itype = 3;
3039 x = 1+(rkappa-BKMXX1)*FBKX1;
3054 AC[1] = U1[1] + U1[2]*
x + U1[3]*
x2 + U1[5]*
y + U1[6]*y2 + U1[7]*y3 +
3055 U1[8]*
xy + U1[10]*p3 + U1[12]*q3 + U1[13]*pq;
3056 AC[2] = U2[1] + U2[2]*
x + U2[3]*
x2 + U2[5]*
y + U2[6]*y2 + U2[7]*y3 +
3057 U2[8]*
xy + U2[9]*p2 + U2[10]*p3 + U2[12]*q3 + U2[13]*pq;
3058 AC[3] = U3[1] + U3[2]*
x + U3[3]*
x2 + U3[5]*
y + U3[6]*y2 + U3[7]*y3 +
3059 U3[8]*
xy + U3[9]*p2 + U3[10]*p3 + U3[11]*q2 + U3[12]*q3 + U3[13]*pq;
3060 AC[4] = U4[1] + U4[2]*
x + U4[3]*
x2 + U4[4]*
x3 + U4[5]*
y + U4[6]*y2 + U4[7]*y3 +
3061 U4[8]*
xy + U4[9]*p2 + U4[10]*p3 + U4[11]*q2 + U4[12]*q3;
3062 AC[5] = U5[1] + U5[2]*
x + U5[3]*
x2 + U5[4]*
x3 + U5[5]*
y + U5[6]*y2 + U5[7]*y3 +
3063 U5[8]*
xy + U5[10]*p3 + U5[11]*q2 + U5[12]*q3 + U5[13]*pq;
3064 AC[6] = U6[1] + U6[2]*
x + U6[3]*
x2 + U6[4]*
x3 + U6[5]*
y + U6[7]*y3 +
3065 U6[8]*
xy + U6[9]*p2 + U6[10]*p3 + U6[12]*q3 + U6[13]*pq;
3066 AC[7] = U7[1] + U7[2]*
x + U7[3]*
x2 + U7[4]*
x3 + U7[5]*
y + U7[6]*y2 + U7[8]*
xy;
3068 AC[8] = U8[1] + U8[2]*
x + U8[3]*
x2 + U8[4]*
x3 + U8[5]*
y + U8[6]*y2 + U8[7]*y3 +
3069 U8[8]*
xy + U8[9]*p2 + U8[10]*p3 + U8[11]*q2 + U8[13]*pq;
3073 AC[9] = (AC[8] - AC[0])/npt;
3076 x = (AC[7]-AC[8])/(AC[7]*AC[8]);
3079 AC[11] = p2*(AC[1]*
TMath::Exp(-AC[2]*(AC[7]+AC[5]*p2)-
3080 AC[3]*
TMath::Exp(-AC[4]*(AC[7]+AC[6]*p2)))-0.045*
y/AC[7])/(1+
x*
y*AC[7]);
3081 AC[12] = (0.045+
x*AC[11])*
y;
3083 if (itype == 4) AC[13] = 0.995/
LandauI(AC[8]);
3085 if (mode==0)
return;
3092 for (k=1; k<=npt; k++) {
3095 WCM[k] = WCM[k-1] + fl + fu;
3099 for (k=1; k<=npt; k++)
3109 if (rlam < AC[0] || rlam > AC[8])
3116 x = (rlam +
HC[0])*
HC[1];
3119 for (k=2; k<=8; k++) {
3121 h[k+1] =
x*
h[k]-fn*
h[k-1];
3124 for (k=2; k<=6; k++)
3128 else if (itype == 2) {
3132 else if (itype == 3) {
3138 v = (AC[11]*
x + AC[12])*
x;
3141 else if (itype == 4) {
3149 #ifdef R__HAS_VECCORE