Garfield++ 3.0
A toolkit for the detailed simulation of particle detectors based on ionisation measurement in gases and semiconductors
Loading...
Searching...
No Matches
Namespaces | Functions
Garfield::Numerics Namespace Reference

Collection of numerical routines. More...

Namespaces

namespace  CERNLIB
 Linear algebra routines from CERNLIB.
 
namespace  QUADPACK
 

Functions

double Legendre (const unsigned int n, const double x)
 Legendre polynomials.
 
double BesselI0S (const double xx)
 
double BesselI1S (const double xx)
 
double BesselK0S (const double xx)
 
double BesselK0L (const double xx)
 
double BesselK1S (const double xx)
 
double BesselK1L (const double xx)
 
double Divdif (const std::vector< double > &f, const std::vector< double > &a, int nn, double x, int mm)
 
bool Boxin2 (const std::vector< std::vector< double > > &value, const std::vector< double > &xAxis, const std::vector< double > &yAxis, const int nx, const int ny, const double xx, const double yy, double &f, const int iOrder)
 
bool Boxin3 (const std::vector< std::vector< std::vector< double > > > &value, const std::vector< double > &xAxis, const std::vector< double > &yAxis, const std::vector< double > &zAxis, const int nx, const int ny, const int nz, const double xx, const double yy, const double zz, double &f, const int iOrder)
 
bool LeastSquaresFit (std::function< double(double, const std::vector< double > &)> f, std::vector< double > &par, std::vector< double > &epar, const std::vector< double > &x, const std::vector< double > &y, const std::vector< double > &ey, const unsigned int nMaxIter, const double diff, double &chi2, const double eps, const bool debug, const bool verbose)
 Least-squares minimisation.
 

Detailed Description

Collection of numerical routines.

Function Documentation

◆ BesselI0S()

double Garfield::Numerics::BesselI0S ( const double  xx)
inline

Modified Bessel functions. Series expansions from Abramowitz and Stegun.

Definition at line 135 of file Numerics.hh.

135 {
136 const double y = xx / 3.75;
137 const double y2 = y * y;
138 return 1. + 3.5156229 * y2 + 3.0899424 * y2 * y2 +
139 1.2067492 * pow(y2, 3) + 0.2659732 * pow(y2, 4) +
140 0.0360768 * pow(y2, 5) + 0.0045813 * pow(y2, 6);
141}

Referenced by BesselK0S().

◆ BesselI1S()

double Garfield::Numerics::BesselI1S ( const double  xx)
inline

Definition at line 143 of file Numerics.hh.

143 {
144 const double y = xx / 3.75;
145 const double y2 = y * y;
146 return xx *
147 (0.5 + 0.87890594 * y2 + 0.51498869 * y2 * y2 +
148 0.15084934 * pow(y2, 3) + 0.02658733 * pow(y2, 4) +
149 0.00301532 * pow(y2, 5) + 0.00032411 * pow(y2, 6));
150}

Referenced by BesselK1S().

◆ BesselK0L()

double Garfield::Numerics::BesselK0L ( const double  xx)
inline

Definition at line 161 of file Numerics.hh.

161 {
162 const double y = 2. / xx;
163 return (exp(-xx) / sqrt(xx)) *
164 (1.25331414 - 0.07832358 * y + 0.02189568 * y * y -
165 0.01062446 * pow(y, 3) + 0.00587872 * pow(y, 4) -
166 0.00251540 * pow(y, 5) + 0.00053208 * pow(y, 6));
167}

◆ BesselK0S()

double Garfield::Numerics::BesselK0S ( const double  xx)
inline

Definition at line 152 of file Numerics.hh.

152 {
153 const double y = xx / 2.;
154 const double y2 = y * y;
155 return -log(y) * BesselI0S(xx) - 0.57721566 +
156 0.42278420 * y2 + 0.23069756 * y2 * y2 +
157 0.03488590 * pow(y2, 3) + 0.00262698 * pow(y2, 4) +
158 0.00010750 * pow(y2, 5) + 0.00000740 * pow(y2, 6);
159}
Garfield::Numerics::BesselI0S
double BesselI0S(const double xx)
Definition: Numerics.hh:135

◆ BesselK1L()

double Garfield::Numerics::BesselK1L ( const double  xx)
inline

Definition at line 179 of file Numerics.hh.

179 {
180 const double y = 2. / xx;
181 return (exp(-xx) / sqrt(xx)) *
182 (1.25331414 + 0.23498619 * y - 0.03655620 * y * y +
183 0.01504268 * pow(y, 3) - 0.00780353 * pow(y, 4) +
184 0.00325614 * pow(y, 5) - 0.00068245 * pow(y, 6));
185}

◆ BesselK1S()

double Garfield::Numerics::BesselK1S ( const double  xx)
inline

Definition at line 169 of file Numerics.hh.

169 {
170 const double y = xx / 2.;
171 const double y2 = y * y;
172 return log(y) * BesselI1S(xx) +
173 (1. / xx) *
174 (1. + 0.15443144 * y2 - 0.67278579 * y2 * y2 -
175 0.18156897 * pow(y2, 3) - 0.01919402 * pow(y2, 4) -
176 0.00110404 * pow(y2, 5) - 0.00004686 * pow(y2, 6));
177}
Garfield::Numerics::BesselI1S
double BesselI1S(const double xx)
Definition: Numerics.hh:143

◆ Boxin2()

bool Garfield::Numerics::Boxin2 ( const std::vector< std::vector< double > > &  value,
const std::vector< double > &  xAxis,
const std::vector< double > &  yAxis,
const int  nx,
const int  ny,
const double  xx,
const double  yy,
double &  f,
const int  iOrder 
)

Interpolation of order 1 and 2 in an irregular rectangular two-dimensional grid.

Definition at line 1304 of file Numerics.cc.

1307 {
1308 //-----------------------------------------------------------------------
1309 // BOXIN2 - Interpolation of order 1 and 2 in an irregular rectangular
1310 // 2-dimensional grid.
1311 //-----------------------------------------------------------------------
1312 int iX0 = 0, iX1 = 0;
1313 int iY0 = 0, iY1 = 0;
1314 std::array<double, 3> fX;
1315 std::array<double, 3> fY;
1316 f = 0.;
1317 // Ensure we are in the grid.
1318 if ((xAxis[nx - 1] - x) * (x - xAxis[0]) < 0 ||
1319 (yAxis[ny - 1] - y) * (y - yAxis[0]) < 0) {
1320 std::cerr << "Boxin2: Point not in the grid; no interpolation.\n";
1321 return false;
1322 }
1323 // Make sure we have enough points.
1324 if (iOrder < 0 || iOrder > 2) {
1325 std::cerr << "Boxin2: Incorrect order; no interpolation.\n";
1326 return false;
1327 } else if (nx < 1 || ny < 1) {
1328 std::cerr << "Boxin2: Incorrect number of points; no interpolation.\n";
1329 return false;
1330 }
1331 if (iOrder == 0 || nx <= 1) {
1332 // Zeroth order interpolation in x.
1333 // Find the nearest node.
1334 double dist = fabs(x - xAxis[0]);
1335 int iNode = 0;
1336 for (int i = 1; i < nx; i++) {
1337 if (fabs(x - xAxis[i]) < dist) {
1338 dist = fabs(x - xAxis[i]);
1339 iNode = i;
1340 }
1341 }
1342 // Set the summing range.
1343 iX0 = iNode;
1344 iX1 = iNode;
1345 // Establish the shape functions.
1346 fX = {1., 0., 0.};
1347 } else if (iOrder == 1 || nx <= 2) {
1348 // First order interpolation in x.
1349 // Find the grid segment containing this point.
1350 int iGrid = 0;
1351 for (int i = 1; i < nx; i++) {
1352 if ((xAxis[i - 1] - x) * (x - xAxis[i]) >= 0.) {
1353 iGrid = i;
1354 }
1355 }
1356 const double x0 = xAxis[iGrid - 1];
1357 const double x1 = xAxis[iGrid];
1358 // Ensure there won't be divisions by zero.
1359 if (x1 == x0) {
1360 std::cerr << "Boxin2: Incorrect grid; no interpolation.\n";
1361 return false;
1362 }
1363 // Compute local coordinates.
1364 const double xL = (x - x0) / (x1 - x0);
1365 // Set the summing range.
1366 iX0 = iGrid - 1;
1367 iX1 = iGrid;
1368 // Set the shape functions.
1369 fX = {1. - xL, xL, 0.};
1370 } else if (iOrder == 2) {
1371 // Second order interpolation in x.
1372 // Find the nearest node and the grid segment.
1373 double dist = fabs(x - xAxis[0]);
1374 int iNode = 0;
1375 for (int i = 1; i < nx; ++i) {
1376 if (fabs(x - xAxis[i]) < dist) {
1377 dist = fabs(x - xAxis[i]);
1378 iNode = i;
1379 }
1380 }
1381 // Find the nearest fitting 2x2 matrix.
1382 int iGrid = std::max(1, std::min(nx - 2, iNode));
1383 const double x0 = xAxis[iGrid - 1];
1384 const double x1 = xAxis[iGrid];
1385 const double x2 = xAxis[iGrid + 1];
1386 // Ensure there won't be divisions by zero.
1387 if (x2 == x0) {
1388 std::cerr << "Boxin2: Incorrect grid; no interpolation.\n";
1389 return false;
1390 }
1391 // Compute the alpha and local coordinate for this grid segment.
1392 const double xAlpha = (x1 - x0) / (x2 - x0);
1393 const double xL = (x - x0) / (x2 - x0);
1394 // Ensure there won't be divisions by zero.
1395 if (xAlpha <= 0 || xAlpha >= 1) {
1396 std::cerr << "Boxin2: Incorrect grid; no interpolation.\n";
1397 return false;
1398 }
1399 // Set the summing range.
1400 iX0 = iGrid - 1;
1401 iX1 = iGrid + 1;
1402 // Set the shape functions.
1403 const double xL2 = xL * xL;
1404 fX[0] = xL2 / xAlpha - xL * (1. + xAlpha) / xAlpha + 1.;
1405 fX[1] = (xL2 - xL) / (xAlpha * xAlpha - xAlpha);
1406 fX[2] = (xL2 - xL * xAlpha) / (1. - xAlpha);
1407 }
1408 if (iOrder == 0 || ny <= 1) {
1409 // Zeroth order interpolation in y.
1410 // Find the nearest node.
1411 double dist = fabs(y - yAxis[0]);
1412 int iNode = 0;
1413 for (int i = 1; i < ny; i++) {
1414 if (fabs(y - yAxis[i]) < dist) {
1415 dist = fabs(y - yAxis[i]);
1416 iNode = i;
1417 }
1418 }
1419 // Set the summing range.
1420 iY0 = iNode;
1421 iY1 = iNode;
1422 // Establish the shape functions.
1423 fY = {1., 0., 0.};
1424 } else if (iOrder == 1 || ny <= 2) {
1425 // First order interpolation in y.
1426 // Find the grid segment containing this point.
1427 int iGrid = 0;
1428 for (int i = 1; i < ny; ++i) {
1429 if ((yAxis[i - 1] - y) * (y - yAxis[i]) >= 0) {
1430 iGrid = i;
1431 }
1432 }
1433 const double y0 = yAxis[iGrid - 1];
1434 const double y1 = yAxis[iGrid];
1435 // Ensure there won't be divisions by zero.
1436 if (y1 == y0) {
1437 std::cerr << "Boxin2: Incorrect grid; no interpolation.\n";
1438 return false;
1439 }
1440 // Compute local coordinates.
1441 const double yL = (y - y0) / (y1 - y0);
1442 // Set the summing range.
1443 iY0 = iGrid - 1;
1444 iY1 = iGrid;
1445 // Set the shape functions.
1446 fY = {1. - yL, yL, 0.};
1447 } else if (iOrder == 2) {
1448 // Second order interpolation in y.
1449 // Find the nearest node and the grid segment.
1450 double dist = fabs(y - yAxis[0]);
1451 int iNode = 0;
1452 for (int i = 1; i < ny; ++i) {
1453 if (fabs(y - yAxis[i]) < dist) {
1454 dist = fabs(y - yAxis[i]);
1455 iNode = i;
1456 }
1457 }
1458 // Find the nearest fitting 2x2 matrix.
1459 int iGrid = std::max(1, std::min(ny - 2, iNode));
1460 const double y0 = yAxis[iGrid - 1];
1461 const double y1 = yAxis[iGrid];
1462 const double y2 = yAxis[iGrid + 1];
1463 // Ensure there won't be divisions by zero.
1464 if (y2 == y0) {
1465 std::cerr << "Boxin2: Incorrect grid; no interpolation.\n";
1466 return false;
1467 }
1468 // Compute the alpha and local coordinate for this grid segment.
1469 const double yAlpha = (y1 - y0) / (y2 - y0);
1470 const double yL = (y - y0) / (y2 - y0);
1471 // Ensure there won't be divisions by zero.
1472 if (yAlpha <= 0 || yAlpha >= 1) {
1473 std::cerr << "Boxin2: Incorrect grid; no interpolation.\n";
1474 return false;
1475 }
1476 // Set the summing range.
1477 iY0 = iGrid - 1;
1478 iY1 = iGrid + 1;
1479 // Set the shape functions.
1480 const double yL2 = yL * yL;
1481 fY[0] = yL2 / yAlpha - yL * (1. + yAlpha) / yAlpha + 1.;
1482 fY[1] = (yL2 - yL) / (yAlpha * yAlpha - yAlpha);
1483 fY[2] = (yL2 - yL * yAlpha) / (1. - yAlpha);
1484 }
1485
1486 // Sum the shape functions.
1487 for (int i = iX0; i <= iX1; ++i) {
1488 for (int j = iY0; j <= iY1; ++j) {
1489 f += value[i][j] * fX[i - iX0] * fY[j - iY0];
1490 }
1491 }
1492 return true;
1493}
Heed::fabs
DoubleAc fabs(const DoubleAc &f)
Definition: DoubleAc.h:615

◆ Boxin3()

bool Garfield::Numerics::Boxin3 ( const std::vector< std::vector< std::vector< double > > > &  value,
const std::vector< double > &  xAxis,
const std::vector< double > &  yAxis,
const std::vector< double > &  zAxis,
const int  nx,
const int  ny,
const int  nz,
const double  xx,
const double  yy,
const double  zz,
double &  f,
const int  iOrder 
)

Interpolation of order 1 and 2 in an irregular rectangular three-dimensional grid.

Definition at line 1495 of file Numerics.cc.

1499 {
1500 //-----------------------------------------------------------------------
1501 // BOXIN3 - interpolation of order 1 and 2 in an irregular rectangular
1502 // 3-dimensional grid.
1503 //-----------------------------------------------------------------------
1504 int iX0 = 0, iX1 = 0;
1505 int iY0 = 0, iY1 = 0;
1506 int iZ0 = 0, iZ1 = 0;
1507 std::array<double, 4> fX;
1508 std::array<double, 4> fY;
1509 std::array<double, 4> fZ;
1510
1511 f = 0.;
1512 // Ensure we are in the grid.
1513 const double x = std::min(std::max(xx, std::min(xAxis[0], xAxis[nx - 1])),
1514 std::max(xAxis[0], xAxis[nx - 1]));
1515 const double y = std::min(std::max(yy, std::min(yAxis[0], yAxis[ny - 1])),
1516 std::max(yAxis[0], yAxis[ny - 1]));
1517 const double z = std::min(std::max(zz, std::min(zAxis[0], zAxis[nz - 1])),
1518 std::max(zAxis[0], zAxis[nz - 1]));
1519
1520 // Make sure we have enough points.
1521 if (iOrder < 0 || iOrder > 2) {
1522 std::cerr << "Boxin3: Incorrect order; no interpolation.\n";
1523 return false;
1524 } else if (nx < 1 || ny < 1 || nz < 1) {
1525 std::cerr << "Boxin3: Incorrect number of points; no interpolation.\n";
1526 return false;
1527 }
1528 if (iOrder == 0 || nx == 1) {
1529 // Zeroth order interpolation in x.
1530 // Find the nearest node.
1531 double dist = fabs(x - xAxis[0]);
1532 int iNode = 0;
1533 for (int i = 1; i < nx; i++) {
1534 if (fabs(x - xAxis[i]) < dist) {
1535 dist = fabs(x - xAxis[i]);
1536 iNode = i;
1537 }
1538 }
1539 // Set the summing range.
1540 iX0 = iNode;
1541 iX1 = iNode;
1542 // Establish the shape functions.
1543 fX = {1., 0., 0., 0.};
1544 } else if (iOrder == 1 || nx == 2) {
1545 // First order interpolation in x.
1546 // Find the grid segment containing this point.
1547 int iGrid = 0;
1548 for (int i = 1; i < nx; i++) {
1549 if ((xAxis[i - 1] - x) * (x - xAxis[i]) >= 0.) {
1550 iGrid = i;
1551 }
1552 }
1553 const double x0 = xAxis[iGrid - 1];
1554 const double x1 = xAxis[iGrid];
1555 // Ensure there won't be divisions by zero.
1556 if (x1 == x0) {
1557 std::cerr << "Boxin3: Incorrect grid; no interpolation.\n";
1558 return false;
1559 }
1560 // Compute local coordinates.
1561 const double xL = (x - x0) / (x1 - x0);
1562 // Set the summing range.
1563 iX0 = iGrid - 1;
1564 iX1 = iGrid;
1565 // Set the shape functions.
1566 fX = {1. - xL, xL, 0., 0.};
1567 } else if (iOrder == 2) {
1568 // Second order interpolation in x.
1569 // Find the grid segment containing this point.
1570 int iGrid = 0;
1571 for (int i = 1; i < nx; i++) {
1572 if ((xAxis[i - 1] - x) * (x - xAxis[i]) >= 0.) {
1573 iGrid = i;
1574 }
1575 }
1576 if (iGrid == 1) {
1577 iX0 = iGrid - 1;
1578 iX1 = iGrid + 1;
1579 const double x0 = xAxis[iX0];
1580 const double x1 = xAxis[iX0 + 1];
1581 const double x2 = xAxis[iX0 + 2];
1582 if (x0 == x1 || x0 == x2 || x1 == x2) {
1583 std::cerr << "Boxin3: One or more grid points in x coincide.\n"
1584 << " No interpolation.\n";
1585 return false;
1586 }
1587 fX[0] = (x - x1) * (x - x2) / ((x0 - x1) * (x0 - x2));
1588 fX[1] = (x - x0) * (x - x2) / ((x1 - x0) * (x1 - x2));
1589 fX[2] = (x - x0) * (x - x1) / ((x2 - x0) * (x2 - x1));
1590 } else if (iGrid == nx - 1) {
1591 iX0 = iGrid - 2;
1592 iX1 = iGrid;
1593 const double x0 = xAxis[iX0];
1594 const double x1 = xAxis[iX0 + 1];
1595 const double x2 = xAxis[iX0 + 2];
1596 if (x0 == x1 || x0 == x2 || x1 == x2) {
1597 std::cerr << "Boxin3: One or more grid points in x coincide.\n"
1598 << " No interpolation.\n";
1599 return false;
1600 }
1601 fX[0] = (x - x1) * (x - x2) / ((x0 - x1) * (x0 - x2));
1602 fX[1] = (x - x0) * (x - x2) / ((x1 - x0) * (x1 - x2));
1603 fX[2] = (x - x0) * (x - x1) / ((x2 - x0) * (x2 - x1));
1604 } else {
1605 iX0 = iGrid - 2;
1606 iX1 = iGrid + 1;
1607 const double x0 = xAxis[iX0];
1608 const double x1 = xAxis[iX0 + 1];
1609 const double x2 = xAxis[iX0 + 2];
1610 const double x3 = xAxis[iX0 + 3];
1611 if (x0 == x1 || x0 == x2 || x0 == x3 ||
1612 x1 == x2 || x1 == x3 || x2 == x3) {
1613 std::cerr << "Boxin3: One or more grid points in x coincide.\n"
1614 << " No interpolation.\n";
1615 return false;
1616 }
1617 // Compute the local coordinate for this grid segment.
1618 const double xL = (x - x1) / (x2 - x1);
1619 fX[0] = ((x - x1) * (x - x2) / ((x0 - x1) * (x0 - x2))) * (1. - xL);
1620 fX[1] = ((x - x0) * (x - x2) / ((x1 - x0) * (x1 - x2))) * (1. - xL) +
1621 ((x - x2) * (x - x3) / ((x1 - x2) * (x1 - x3))) * xL;
1622 fX[2] = ((x - x0) * (x - x1) / ((x2 - x0) * (x2 - x1))) * (1. - xL) +
1623 ((x - x1) * (x - x3) / ((x2 - x1) * (x2 - x3))) * xL;
1624 fX[3] = ((x - x1) * (x - x2) / ((x3 - x1) * (x3 - x2))) * xL;
1625 }
1626 }
1627
1628 if (iOrder == 0 || ny == 1) {
1629 // Zeroth order interpolation in y.
1630 // Find the nearest node.
1631 double dist = fabs(y - yAxis[0]);
1632 int iNode = 0;
1633 for (int i = 1; i < ny; i++) {
1634 if (fabs(y - yAxis[i]) < dist) {
1635 dist = fabs(y - yAxis[i]);
1636 iNode = i;
1637 }
1638 }
1639 // Set the summing range.
1640 iY0 = iNode;
1641 iY1 = iNode;
1642 // Establish the shape functions.
1643 fY = {1., 0., 0., 0.};
1644 } else if (iOrder == 1 || ny == 2) {
1645 // First order interpolation in y.
1646 // Find the grid segment containing this point.
1647 int iGrid = 0;
1648 for (int i = 1; i < ny; i++) {
1649 if ((yAxis[i - 1] - y) * (y - yAxis[i]) >= 0.) {
1650 iGrid = i;
1651 }
1652 }
1653 // Ensure there won't be divisions by zero.
1654 const double y0 = yAxis[iGrid - 1];
1655 const double y1 = yAxis[iGrid];
1656 if (y1 == y0) {
1657 std::cerr << "Boxin3: Incorrect grid; no interpolation.\n";
1658 return false;
1659 }
1660 // Compute local coordinates.
1661 const double yL = (y - y0) / (y1 - y0);
1662 // Set the summing range.
1663 iY0 = iGrid - 1;
1664 iY1 = iGrid;
1665 // Set the shape functions.
1666 fY = {1. - yL, yL, 0., 0.};
1667 } else if (iOrder == 2) {
1668 // Second order interpolation in y.
1669 // Find the grid segment containing this point.
1670 int iGrid = 0;
1671 for (int i = 1; i < ny; i++) {
1672 if ((yAxis[i - 1] - y) * (y - yAxis[i]) >= 0.) {
1673 iGrid = i;
1674 }
1675 }
1676 if (iGrid == 1) {
1677 iY0 = iGrid - 1;
1678 iY1 = iGrid + 1;
1679 const double y0 = yAxis[iY0];
1680 const double y1 = yAxis[iY0 + 1];
1681 const double y2 = yAxis[iY0 + 2];
1682 if (y0 == y1 || y0 == y2 || y1 == y2) {
1683 std::cerr << "Boxin3: One or more grid points in y coincide.\n"
1684 << " No interpolation.\n";
1685 return false;
1686 }
1687 fY[0] = (y - y1) * (y - y2) / ((y0 - y1) * (y0 - y2));
1688 fY[1] = (y - y0) * (y - y2) / ((y1 - y0) * (y1 - y2));
1689 fY[2] = (y - y0) * (y - y1) / ((y2 - y0) * (y2 - y1));
1690 } else if (iGrid == ny - 1) {
1691 iY0 = iGrid - 2;
1692 iY1 = iGrid;
1693 const double y0 = yAxis[iY0];
1694 const double y1 = yAxis[iY0 + 1];
1695 const double y2 = yAxis[iY0 + 2];
1696 if (y0 == y1 || y0 == y2 || y1 == y2) {
1697 std::cerr << "Boxin3: One or more grid points in y coincide.\n"
1698 << " No interpolation.\n";
1699 return false;
1700 }
1701 fY[0] = (y - y1) * (y - y2) / ((y0 - y1) * (y0 - y2));
1702 fY[1] = (y - y0) * (y - y2) / ((y1 - y0) * (y1 - y2));
1703 fY[2] = (y - y0) * (y - y1) / ((y2 - y0) * (y2 - y1));
1704 } else {
1705 iY0 = iGrid - 2;
1706 iY1 = iGrid + 1;
1707 const double y0 = yAxis[iY0];
1708 const double y1 = yAxis[iY0 + 1];
1709 const double y2 = yAxis[iY0 + 2];
1710 const double y3 = yAxis[iY0 + 3];
1711 if (y0 == y1 || y0 == y2 || y0 == y3 ||
1712 y1 == y2 || y1 == y3 || y2 == y3) {
1713 std::cerr << "Boxin3: One or more grid points in y coincide.\n"
1714 << " No interpolation.\n";
1715 return false;
1716 }
1717 // Compute the local coordinate for this grid segment.
1718 const double yL = (y - y1) / (y2 - y1);
1719 fY[0] = ((y - y1) * (y - y2) / ((y0 - y1) * (y0 - y2))) * (1. - yL);
1720 fY[1] = ((y - y0) * (y - y2) / ((y1 - y0) * (y1 - y2))) * (1. - yL) +
1721 ((y - y2) * (y - y3) / ((y1 - y2) * (y1 - y3))) * yL;
1722 fY[2] = ((y - y0) * (y - y1) / ((y2 - y0) * (y2 - y1))) * (1. - yL) +
1723 ((y - y1) * (y - y3) / ((y2 - y1) * (y2 - y3))) * yL;
1724 fY[3] = ((y - y1) * (y - y2) / ((y3 - y1) * (y3 - y2))) * yL;
1725 }
1726 }
1727
1728 if (iOrder == 0 || nz == 1) {
1729 // Zeroth order interpolation in z.
1730 // Find the nearest node.
1731 double dist = fabs(z - zAxis[0]);
1732 int iNode = 0;
1733 for (int i = 1; i < nz; i++) {
1734 if (fabs(z - zAxis[i]) < dist) {
1735 dist = fabs(z - zAxis[i]);
1736 iNode = i;
1737 }
1738 }
1739 // Set the summing range.
1740 iZ0 = iNode;
1741 iZ1 = iNode;
1742 // Establish the shape functions.
1743 fZ = {1., 0., 0., 0.};
1744 } else if (iOrder == 1 || nz == 2) {
1745 // First order interpolation in z.
1746 // Find the grid segment containing this point.
1747 int iGrid = 0;
1748 for (int i = 1; i < nz; i++) {
1749 if ((zAxis[i - 1] - z) * (z - zAxis[i]) >= 0.) {
1750 iGrid = i;
1751 }
1752 }
1753 const double z0 = zAxis[iGrid - 1];
1754 const double z1 = zAxis[iGrid];
1755 // Ensure there won't be divisions by zero.
1756 if (z1 == z0) {
1757 std::cerr << "Boxin3: Incorrect grid; no interpolation.\n";
1758 return false;
1759 }
1760 // Compute local coordinates.
1761 const double zL = (z - z0) / (z1 - z0);
1762 // Set the summing range.
1763 iZ0 = iGrid - 1;
1764 iZ1 = iGrid;
1765 // Set the shape functions.
1766 fZ = {1. - zL, zL, 0., 0.};
1767 } else if (iOrder == 2) {
1768 // Second order interpolation in z.
1769 // Find the grid segment containing this point.
1770 int iGrid = 0;
1771 for (int i = 1; i < nz; i++) {
1772 if ((zAxis[i - 1] - z) * (z - zAxis[i]) >= 0.) {
1773 iGrid = i;
1774 }
1775 }
1776 if (iGrid == 1) {
1777 iZ0 = iGrid - 1;
1778 iZ1 = iGrid + 1;
1779 const double z0 = zAxis[iZ0];
1780 const double z1 = zAxis[iZ0 + 1];
1781 const double z2 = zAxis[iZ0 + 2];
1782 if (z0 == z1 || z0 == z2 || z1 == z2) {
1783 std::cerr << "Boxin3: One or more grid points in z coincide.\n"
1784 << " No interpolation.\n";
1785 return false;
1786 }
1787 fZ[0] = (z - z1) * (z - z2) / ((z0 - z1) * (z0 - z2));
1788 fZ[1] = (z - z0) * (z - z2) / ((z1 - z0) * (z1 - z2));
1789 fZ[2] = (z - z0) * (z - z1) / ((z2 - z0) * (z2 - z1));
1790 } else if (iGrid == nz - 1) {
1791 iZ0 = iGrid - 2;
1792 iZ1 = iGrid;
1793 const double z0 = zAxis[iZ0];
1794 const double z1 = zAxis[iZ0 + 1];
1795 const double z2 = zAxis[iZ0 + 2];
1796 if (z0 == z1 || z0 == z2 || z1 == z2) {
1797 std::cerr << "Boxin3: One or more grid points in z coincide.\n"
1798 << " No interpolation.\n";
1799 return false;
1800 }
1801 fZ[0] = (z - z1) * (z - z2) / ((z0 - z1) * (z0 - z2));
1802 fZ[1] = (z - z0) * (z - z2) / ((z1 - z0) * (z1 - z2));
1803 fZ[2] = (z - z0) * (z - z1) / ((z2 - z0) * (z2 - z1));
1804 } else {
1805 iZ0 = iGrid - 2;
1806 iZ1 = iGrid + 1;
1807 const double z0 = zAxis[iZ0];
1808 const double z1 = zAxis[iZ0 + 1];
1809 const double z2 = zAxis[iZ0 + 2];
1810 const double z3 = zAxis[iZ0 + 3];
1811 if (z0 == z1 || z0 == z2 || z0 == z3 ||
1812 z1 == z2 || z1 == z3 || z2 == z3) {
1813 std::cerr << "Boxin3: One or more grid points in z coincide.\n"
1814 << " No interpolation.\n";
1815 return false;
1816 }
1817 // Compute the local coordinate for this grid segment.
1818 const double zL = (z - z1) / (z2 - z1);
1819 fZ[0] = ((z - z1) * (z - z2) / ((z0 - z1) * (z0 - z2))) * (1. - zL);
1820 fZ[1] = ((z - z0) * (z - z2) / ((z1 - z0) * (z1 - z2))) * (1. - zL) +
1821 ((z - z2) * (z - z3) / ((z1 - z2) * (z1 - z3))) * zL;
1822 fZ[2] = ((z - z0) * (z - z1) / ((z2 - z0) * (z2 - z1))) * (1. - zL) +
1823 ((z - z1) * (z - z3) / ((z2 - z1) * (z2 - z3))) * zL;
1824 fZ[3] = ((z - z1) * (z - z2) / ((z3 - z1) * (z3 - z2))) * zL;
1825 }
1826 }
1827
1828 for (int i = iX0; i <= iX1; ++i) {
1829 for (int j = iY0; j <= iY1; ++j) {
1830 for (int k = iZ0; k <= iZ1; ++k) {
1831 f += value[i][j][k] * fX[i - iX0] * fY[j - iY0] * fZ[k - iZ0];
1832 }
1833 }
1834 }
1835 return true;
1836}

Referenced by Garfield::Medium::Interpolate().

◆ Divdif()

double Garfield::Numerics::Divdif ( const std::vector< double > &  f,
const std::vector< double > &  a,
int  nn,
double  x,
int  mm 
)

C++ version of DIVDIF (CERN program library E105) which performs tabular interpolation using symmetrically placed argument points.

Definition at line 1206 of file Numerics.cc.

1207 {
1208
1209 double t[20], d[20];
1210
1211 // Check the arguments.
1212 if (nn < 2) {
1213 std::cerr << "Divdif: Array length < 2.\n";
1214 return 0.;
1215 }
1216 if (mm < 1) {
1217 std::cerr << "Divdif: Interpolation order < 1.\n";
1218 return 0.;
1219 }
1220
1221 // Deal with the case that X is located at first or last point.
1222 const double tol = 1.e-6 * (fabs(a[0]) + fabs(a[nn - 1]));
1223 if (fabs(x - a[0]) < tol) return f[0];
1224 if (fabs(x - a[nn - 1]) < tol) return f[nn - 1];
1225
1226 // Find subscript IX of X in array A.
1227 constexpr int mmax = 10;
1228 const int m = std::min({mm, mmax, nn - 1});
1229 const int mplus = m + 1;
1230 int ix = 0;
1231 int iy = nn + 1;
1232 if (a[0] > a[nn - 1]) {
1233 // Search decreasing arguments.
1234 do {
1235 const int mid = (ix + iy) / 2;
1236 if (x > a[mid - 1]) {
1237 iy = mid;
1238 } else {
1239 ix = mid;
1240 }
1241 } while (iy - ix > 1);
1242 } else {
1243 // Search increasing arguments.
1244 do {
1245 const int mid = (ix + iy) / 2;
1246 if (x < a[mid - 1]) {
1247 iy = mid;
1248 } else {
1249 ix = mid;
1250 }
1251 } while (iy - ix > 1);
1252 }
1253 // Copy reordered interpolation points into (T[I],D[I]), setting
1254 // EXTRA to True if M+2 points to be used.
1255 int npts = m + 2 - (m % 2);
1256 int ip = 0;
1257 int l = 0;
1258 do {
1259 const int isub = ix + l;
1260 if ((1 > isub) || (isub > nn)) {
1261 // Skip point.
1262 npts = mplus;
1263 } else {
1264 // Insert point.
1265 ip++;
1266 t[ip - 1] = a[isub - 1];
1267 d[ip - 1] = f[isub - 1];
1268 }
1269 if (ip < npts) {
1270 l = -l;
1271 if (l >= 0) ++l;
1272 }
1273 } while (ip < npts);
1274
1275 const bool extra = npts != mplus;
1276 // Replace d by the leading diagonal of a divided-difference table,
1277 // supplemented by an extra line if EXTRA is True.
1278 for (l = 1; l <= m; l++) {
1279 if (extra) {
1280 const int isub = mplus - l;
1281 d[m + 1] = (d[m + 1] - d[m - 1]) / (t[m + 1] - t[isub - 1]);
1282 }
1283 int i = mplus;
1284 for (int j = l; j <= m; j++) {
1285 const int isub = i - l;
1286 d[i - 1] = (d[i - 1] - d[i - 2]) / (t[i - 1] - t[isub - 1]);
1287 i--;
1288 }
1289 }
1290 // Evaluate the Newton interpolation formula at X, averaging two values
1291 // of last difference if EXTRA is True.
1292 double sum = d[mplus - 1];
1293 if (extra) {
1294 sum = 0.5 * (sum + d[m + 1]);
1295 }
1296 int j = m;
1297 for (l = 1; l <= m; l++) {
1298 sum = d[j - 1] + (x - t[j - 1]) * sum;
1299 j--;
1300 }
1301 return sum;
1302}

Referenced by Garfield::Sensor::ComputeThresholdCrossings(), and Garfield::Medium::Interpolate1D().

◆ LeastSquaresFit()

bool Garfield::Numerics::LeastSquaresFit ( std::function< double(double, const std::vector< double > &)>  f,
std::vector< double > &  par,
std::vector< double > &  epar,
const std::vector< double > &  x,
const std::vector< double > &  y,
const std::vector< double > &  ey,
const unsigned int  nMaxIter,
const double  diff,
double &  chi2,
const double  eps,
const bool  debug,
const bool  verbose 
)

Least-squares minimisation.

Definition at line 1838 of file Numerics.cc.

1844 {
1845
1846 //-----------------------------------------------------------------------
1847 // LSQFIT - Subroutine fitting the parameters A in the routine F to
1848 // the data points (X,Y) using a least squares method.
1849 // Translated from an Algol routine written by Geert Jan van
1850 // Oldenborgh and Rob Veenhof, based on Stoer + Bulirsch.
1851 // VARIABLES : F( . ,A,VAL) : Subroutine to be fitted.
1852 // (X,Y) : Input data.
1853 // D : Derivative matrix.
1854 // R : Difference vector between Y and F(X,A).
1855 // S : Correction vector for A.
1856 // EPSDIF : Used for differentiating.
1857 // EPS : Numerical resolution.
1858 // (Last updated on 23/ 5/11.)
1859 //-----------------------------------------------------------------------
1860
1861 const unsigned int n = par.size();
1862 const unsigned int m = x.size();
1863 // Make sure that the # degrees of freedom < the number of data points.
1864 if (n > m) {
1865 std::cerr << "LeastSquaresFit: Number of parameters to be varied\n"
1866 << " exceeds the number of data points.\n";
1867 return false;
1868 }
1869
1870 // Check the errors.
1871 if (*std::min_element(ey.cbegin(), ey.cend()) <= 0.) {
1872 std::cerr << "LeastSquaresFit: Not all errors are > 0; no fit done.\n";
1873 return false;
1874 }
1875 chi2 = 0.;
1876 // Largest difference.
1877 double diffc = -1.;
1878 // Initialise the difference vector R.
1879 std::vector<double> r(m, 0.);
1880 for (unsigned int i = 0; i < m; ++i) {
1881 // Compute initial residuals.
1882 r[i] = (y[i] - f(x[i], par)) / ey[i];
1883 // Compute initial maximum difference.
1884 diffc = std::max(std::abs(r[i]), diffc);
1885 // And compute initial chi2.
1886 chi2 += r[i] * r[i];
1887 }
1888 if (debug) {
1889 std::cout << " Input data points:\n"
1890 << " X Y Y - F(X)\n";
1891 for (unsigned int i = 0; i < m; ++i) {
1892 std::printf(" %9u %15.8e %15.8e %15.8e\n", i, x[i], y[i], r[i]);
1893 }
1894 std::cout << " Initial values of the fit parameters:\n"
1895 << " Parameter Value\n";
1896 for (unsigned int i = 0; i < n; ++i) {
1897 std::printf(" %9u %15.8e\n", i, par[i]);
1898 }
1899 }
1900 if (verbose) {
1901 std::cout << " MINIMISATION SUMMARY\n"
1902 << " Initial situation:\n";
1903 std::printf(" Largest difference between fit and target: %15.8e\n",
1904 diffc);
1905 std::printf(" Sum of squares of these differences: %15.8e\n",
1906 chi2);
1907 }
1908 // Start optimising loop.
1909 bool converged = false;
1910 double chi2L = 0.;
1911 for (unsigned int iter = 1; iter <= nMaxIter; ++iter) {
1912 // Check the stopping criteria: (1) max norm, (2) change in chi-squared.
1913 if ((diffc < diff) || (iter > 1 && std::abs(chi2L - chi2) < eps * chi2)) {
1914 if (debug || verbose) {
1915 if (diffc < diff) {
1916 std::cout << " The maximum difference stopping criterion "
1917 << "is satisfied.\n";
1918 } else {
1919 std::cout << " The relative change in chi-squared has dropped "
1920 << "below the threshold.\n";
1921 }
1922 }
1923 converged = true;
1924 break;
1925 }
1926 // Calculate the derivative matrix.
1927 std::vector<std::vector<double> > d(n, std::vector<double>(m, 0.));
1928 for (unsigned int i = 0; i < n; ++i) {
1929 const double epsdif = eps * (1. + std::abs(par[i]));
1930 par[i] += 0.5 * epsdif;
1931 for (unsigned int j = 0; j < m; ++j) d[i][j] = f(x[j], par);
1932 par[i] -= epsdif;
1933 for (unsigned int j = 0; j < m; ++j) {
1934 d[i][j] = (d[i][j] - f(x[j], par)) / (epsdif * ey[j]);
1935 }
1936 par[i] += 0.5 * epsdif;
1937 }
1938 // Invert the matrix in Householder style.
1939 std::vector<double> colsum(n, 0.);
1940 std::vector<int> pivot(n, 0);
1941 for (unsigned int i = 0; i < n; ++i) {
1942 colsum[i] = std::inner_product(d[i].cbegin(), d[i].cend(),
1943 d[i].cbegin(), 0.);
1944 pivot[i] = i;
1945 }
1946 // Decomposition.
1947 std::vector<double> alpha(n, 0.);
1948 bool singular = false;
1949 for (unsigned int k = 0; k < n; ++k) {
1950 double sigma = colsum[k];
1951 unsigned int jbar = k;
1952 for (unsigned int j = k + 1; j < n; ++j) {
1953 if (sigma < colsum[j]) {
1954 sigma = colsum[j];
1955 jbar = j;
1956 }
1957 }
1958 if (jbar != k) {
1959 // Interchange columns.
1960 std::swap(pivot[k], pivot[jbar]);
1961 std::swap(colsum[k], colsum[jbar]);
1962 std::swap(d[k], d[jbar]);
1963 }
1964 sigma = 0.;
1965 for (unsigned int i = k; i < m; ++i) sigma += d[k][i] * d[k][i];
1966 if (sigma == 0. || sqrt(sigma) < 1.e-8 * std::abs(d[k][k])) {
1967 singular = true;
1968 break;
1969 }
1970 alpha[k] = d[k][k] < 0. ? sqrt(sigma) : -sqrt(sigma);
1971 const double beta = 1. / (sigma - d[k][k] * alpha[k]);
1972 d[k][k] -= alpha[k];
1973 std::vector<double> b(n, 0.);
1974 for (unsigned int j = k + 1; j < n; ++j) {
1975 for (unsigned int i = k; i < n; ++i) b[j] += d[k][i] * d[j][i];
1976 b[j] *= beta;
1977 }
1978 for (unsigned int j = k + 1; j < n; ++j) {
1979 for (unsigned int i = k; i < m; ++i) {
1980 d[j][i] -= d[k][i] * b[j];
1981 colsum[j] -= d[j][k] * d[j][k];
1982 }
1983 }
1984 }
1985 if (singular) {
1986 std::cerr << "LeastSquaresFit: Householder matrix (nearly) singular;\n"
1987 << " no further optimisation.\n"
1988 << " Ensure the function depends on the parameters\n"
1989 << " and try to supply reasonable starting values.\n";
1990 break;
1991 }
1992 // Solve.
1993 for (unsigned int j = 0; j < n; ++j) {
1994 double gamma = 0.;
1995 for (unsigned int i = j; i < m; ++i) gamma += d[j][i] * r[i];
1996 gamma *= 1. / (alpha[j] * d[j][j]);
1997 for (unsigned int i = j; i < m; ++i) r[i] += gamma * d[j][i];
1998 }
1999 std::vector<double> z(n, 0.);
2000 z[n - 1] = r[n - 1] / alpha[n - 1];
2001 for (int i = n - 1; i >= 1; --i) {
2002 double sum = 0.;
2003 for (unsigned int j = i + 1; j <= n; ++j) {
2004 sum += d[j - 1][i - 1] * z[j - 1];
2005 }
2006 z[i - 1] = (r[i - 1] - sum) / alpha[i - 1];
2007 }
2008 // Correction vector.
2009 std::vector<double> s(n, 0.);
2010 for (unsigned int i = 0; i < n; ++i) s[pivot[i]] = z[i];
2011 // Generate some debugging output.
2012 if (debug) {
2013 std::cout << " Correction vector in iteration " << iter << ":\n";
2014 for (unsigned int i = 0; i < n; ++i) {
2015 std::printf(" %5u %15.8e\n", i, s[i]);
2016 }
2017 }
2018 // Add part of the correction vector to the estimate to improve chi2.
2019 chi2L = chi2;
2020 chi2 *= 2;
2021 for (unsigned int i = 0; i < n; ++i) par[i] += s[i] * 2;
2022 for (unsigned int i = 0; i <= 10; ++i) {
2023 if (chi2 <= chi2L) break;
2024 if (std::abs(chi2L - chi2) < eps * chi2) {
2025 if (debug) {
2026 std::cout << " Too little improvement; reduction loop halted.\n";
2027 }
2028 break;
2029 }
2030 chi2 = 0.;
2031 const double scale = 1. / pow(2, i);
2032 for (unsigned int j = 0; j < n; ++j) par[j] -= s[j] * scale;
2033 for (unsigned int j = 0; j < m; ++j) {
2034 r[j] = (y[j] - f(x[j], par)) / ey[j];
2035 chi2 += r[j] * r[j];
2036 }
2037 if (debug) {
2038 std::printf(" Reduction loop %3i: chi2 = %15.8e\n", i, chi2);
2039 }
2040 }
2041 // Calculate the max. norm.
2042 diffc = std::abs(r[0]);
2043 for (unsigned int i = 1; i < m; ++i) {
2044 diffc = std::max(std::abs(r[i]), diffc);
2045 }
2046 // Print some debugging output.
2047 if (debug) {
2048 std::cout << " Values of the fit parameters after iteration "
2049 << iter << "\n Parameter Value\n";
2050 for (unsigned int i = 0; i < n; ++i) {
2051 std::printf(" %9u %15.8e\n", i, par[i]);
2052 }
2053 std::printf(" for which chi2 = %15.8e and diff = %15.8e\n",
2054 chi2, diffc);
2055 } else if (verbose) {
2056 std::printf(" Step %3u: largest deviation = %15.8e, chi2 = %15.8e\n",
2057 iter, diffc, chi2);
2058 }
2059 }
2060 // End of fit, perform error calculation.
2061 if (!converged) {
2062 std::cerr << "LeastSquaresFit: Maximum number of iterations reached.\n";
2063 }
2064 // Calculate the derivative matrix for the final settings.
2065 std::vector<std::vector<double> > d(n, std::vector<double>(m, 0.));
2066 for (unsigned int i = 0; i < n; ++i) {
2067 const double epsdif = eps * (1. + std::abs(par[i]));
2068 par[i] += 0.5 * epsdif;
2069 for (unsigned int j = 0; j < m; ++j) d[i][j] = f(x[j], par);
2070 par[i] -= epsdif;
2071 for (unsigned int j = 0; j < m; ++j) {
2072 d[i][j] = (d[i][j] - f(x[j], par)) / (epsdif * ey[j]);
2073 }
2074 par[i] += 0.5 * epsdif;
2075 }
2076 // Calculate the error matrix.
2077 std::vector<std::vector<double> > cov(n, std::vector<double>(n, 0.));
2078 for (unsigned int i = 0; i < n; ++i) {
2079 for (unsigned int j = 0; j < n; ++j) {
2080 cov[i][j] = std::inner_product(d[i].cbegin(), d[i].cend(),
2081 d[j].cbegin(), 0.);
2082 }
2083 }
2084 // Compute the scaling factor for the errors.
2085 double scale = m > n ? chi2 / (m - n) : 1.;
2086 // Invert it to get the covariance matrix.
2087 epar.assign(n, 0.);
2088 if (Garfield::Numerics::CERNLIB::dinv(n, cov) != 0) {
2089 std::cerr << "LeastSquaresFit: Singular covariance matrix; "
2090 << "no error calculation.\n";
2091 } else {
2092 for (unsigned int i = 0; i < n; ++i) {
2093 for (unsigned int j = 0; j < n; ++j) cov[i][j] *= scale;
2094 epar[i] = sqrt(std::max(0., cov[i][i]));
2095 }
2096 }
2097 // Print results.
2098 if (debug) {
2099 std::cout << " Comparison between input and fit:\n"
2100 << " X Y F(X)\n";
2101 for (unsigned int i = 0; i < m; ++i) {
2102 std::printf(" %5u %15.8e %15.8e %15.8e\n", i, x[i], y[i], f(x[i], par));
2103 }
2104 }
2105 if (verbose) {
2106 std::cout << " Final values of the fit parameters:\n"
2107 << " Parameter Value Error\n";
2108 for (unsigned int i = 0; i < n; ++i) {
2109 std::printf(" %9u %15.8e %15.8e\n", i, par[i], epar[i]);
2110 }
2111 std::cout << " The errors have been scaled by a factor of "
2112 << sqrt(scale) << ".\n";
2113 std::cout << " Covariance matrix:\n";
2114 for (unsigned int i = 0; i < n; ++i) {
2115 for (unsigned int j = 0; j < n; ++j) {
2116 std::printf(" %15.8e", cov[i][j]);
2117 }
2118 std::cout << "\n";
2119 }
2120 std::cout << " Correlation matrix:\n";
2121 for (unsigned int i = 0; i < n; ++i) {
2122 for (unsigned int j = 0; j < n; ++j) {
2123 double cor = 0.;
2124 if (cov[i][i] > 0. && cov[j][j] > 0.) {
2125 cor = cov[i][j] / sqrt(cov[i][i] * cov[j][j]);
2126 }
2127 std::printf(" %15.8e", cor);
2128 }
2129 std::cout << "\n";
2130 }
2131 std::cout << " Minimisation finished.\n";
2132 }
2133 return converged;
2134}
Garfield::Numerics::CERNLIB::dinv
int dinv(const int n, std::vector< std::vector< double > > &a)
Replace square matrix A by its inverse.
Definition: Numerics.cc:787
Heed::pow
DoubleAc pow(const DoubleAc &f, double p)
Definition: DoubleAc.cpp:337
Heed::sqrt
DoubleAc sqrt(const DoubleAc &f)
Definition: DoubleAc.cpp:314

Referenced by Garfield::ComponentAnalyticField::MultipoleMoments().

◆ Legendre()

double Garfield::Numerics::Legendre ( const unsigned int  n,
const double  x 
)
inline

Legendre polynomials.

Definition at line 119 of file Numerics.hh.

119 {
120
121 if (std::abs(x) > 1.) return 0.;
122 double p0 = 1.;
123 double p1 = x;
124 if (n == 0) return p0;
125 if (n == 1) return p1;
126 for (unsigned int k = 1; k < n; ++k) {
127 p0 = ((2 * k + 1) * x * p1 - k * p0) / (k + 1);
128 std::swap(p0, p1);
129 }
130 return p1;
131}

Referenced by Garfield::ComponentAnalyticField::MultipoleMoments().