Garfield++ 5.0
A toolkit for the detailed simulation of particle detectors based on ionisation measurement in gases and semiconductors
Loading...
Searching...
No Matches
Public Types | Public Member Functions | List of all members
Garfield::ComponentAnalyticField Class Reference

#include <ComponentAnalyticField.hh>

+ Inheritance diagram for Garfield::ComponentAnalyticField:

Public Types

enum  Cell {
  A00 , B1X , B1Y , B2X ,
  B2Y , C10 , C2X , C2Y ,
  C30 , D10 , D20 , D30 ,
  D40 , Unknown
}
 

Public Member Functions

 ComponentAnalyticField ()
 Constructor.
 
 ~ComponentAnalyticField ()
 Destructor.
 
void SetMedium (Medium *medium)
 Set the medium inside the cell.
 
void AddWire (const double x, const double y, const double diameter, const double voltage, const std::string &label="", const double length=100., const double tension=50., const double rho=19.3, const int ntrap=5)
 Add a wire at (x, y) .
 
void AddTube (const double radius, const double voltage, const int nEdges, const std::string &label="")
 Add a tube.
 
void AddPlaneX (const double x, const double voltage, const std::string &label="")
 Add a plane at constant x.
 
void AddPlaneY (const double y, const double voltage, const std::string &label="")
 Add a plane at constant y.
 
void AddPlaneR (const double r, const double voltage, const std::string &label="")
 Add a plane at constant radius.
 
void AddPlanePhi (const double phi, const double voltage, const std::string &label="")
 Add a plane at constant phi.
 
void AddStripOnPlaneX (const char direction, const double x, const double smin, const double smax, const std::string &label, const double gap=-1.)
 
void AddStripOnPlaneY (const char direction, const double y, const double smin, const double smax, const std::string &label, const double gap=-1.)
 Add a strip in the x or z direction on an existing plane at constant y.
 
void AddStripOnPlaneR (const char direction, const double r, const double smin, const double smax, const std::string &label, const double gap=-1.)
 Add a strip in the phi or z direction on an existing plane at constant radius.
 
void AddStripOnPlanePhi (const char direction, const double phi, const double smin, const double smax, const std::string &label, const double gap=-1.)
 Add a strip in the r or z direction on an existing plane at constant phi.
 
void AddPixelOnPlaneX (const double x, const double ymin, const double ymax, const double zmin, const double zmax, const std::string &label, const double gap=-1., const double rot=0.)
 
void AddPixelOnPlaneY (const double y, const double xmin, const double xmax, const double zmin, const double zmax, const std::string &label, const double gap=-1., const double rot=0.)
 Add a pixel on an existing plane at constant y.
 
void AddPixelOnPlaneR (const double r, const double phimin, const double phimax, const double zmin, const double zmax, const std::string &label, const double gap=-1.)
 Add a pixel on an existing plane at constant radius.
 
void AddPixelOnPlanePhi (const double phi, const double rmin, const double rmax, const double zmin, const double zmax, const std::string &label, const double gap=-1.)
 Add a pixel on an existing plane at constant phi.
 
void SetPeriodicityX (const double s)
 Set the periodic length [cm] in the x-direction.
 
void SetPeriodicityY (const double s)
 Set the periodic length [cm] in the y-direction.
 
void SetPeriodicityPhi (const double phi)
 Set the periodicity [degree] in phi.
 
bool GetPeriodicityX (double &s)
 Get the periodic length in the x-direction.
 
bool GetPeriodicityY (double &s)
 Get the periodic length in the y-direction.
 
bool GetPeriodicityPhi (double &s)
 Get the periodicity [degree] in phi.
 
void SetCartesianCoordinates ()
 Use Cartesian coordinates (default).
 
void SetPolarCoordinates ()
 Use polar coordinates.
 
bool IsPolar () const
 Are polar coordinates being used?
 
void PrintCell ()
 Print all available information on the cell.
 
void PlotCell (TPad *pad)
 Make a plot of the cell layout.
 
void AddCharge (const double x, const double y, const double z, const double q)
 Add a point charge.
 
void ClearCharges ()
 Remove all point charges.
 
void PrintCharges () const
 Print a list of the point charges.
 
std::string GetCellType ()
 
void AddReadout (const std::string &label, const bool silent=false)
 
void SetNumberOfCellCopies (const unsigned int nfourier)
 
bool MultipoleMoments (const unsigned int iw, const unsigned int order=4, const bool print=false, const bool plot=false, const double rmult=1., const double eps=1.e-4, const unsigned int nMaxIter=20)
 
void EnableDipoleTerms (const bool on=true)
 Request dipole terms be included for each of the wires (default: off).
 
void EnableChargeCheck (const bool on=true)
 Check the quality of the capacitance matrix inversion (default: off).
 
unsigned int GetNumberOfWires () const
 Get the number of wires.
 
bool GetWire (const unsigned int i, double &x, double &y, double &diameter, double &voltage, std::string &label, double &length, double &charge, int &ntrap) const
 Retrieve the parameters of a wire.
 
unsigned int GetNumberOfPlanesX () const
 Get the number of equipotential planes at constant x.
 
unsigned int GetNumberOfPlanesY () const
 Get the number of equipotential planes at constant y.
 
unsigned int GetNumberOfPlanesR () const
 Get the number of equipotential planes at constant radius.
 
unsigned int GetNumberOfPlanesPhi () const
 Get the number of equipotential planes at constant phi.
 
bool GetPlaneX (const unsigned int i, double &x, double &voltage, std::string &label) const
 Retrieve the parameters of a plane at constant x.
 
bool GetPlaneY (const unsigned int i, double &y, double &voltage, std::string &label) const
 Retrieve the parameters of a plane at constant y.
 
bool GetPlaneR (const unsigned int i, double &r, double &voltage, std::string &label) const
 Retrieve the parameters of a plane at constant radius.
 
bool GetPlanePhi (const unsigned int i, double &phi, double &voltage, std::string &label) const
 Retrieve the parameters of a plane at constant phi.
 
bool GetTube (double &r, double &voltage, int &nEdges, std::string &label) const
 Retrieve the tube parameters.
 
bool OptimiseOnTrack (const std::vector< std::string > &groups, const std::string &field_function, const double target, const double x0, const double y0, const double x1, const double y1, const unsigned int nP=20, const bool print=true)
 
bool OptimiseOnGrid (const std::vector< std::string > &groups, const std::string &field_function, const double target, const double x0, const double y0, const double x1, const double y1, const unsigned int nX=10, const unsigned int nY=10, const bool print=true)
 
bool OptimiseOnWires (const std::vector< std::string > &groups, const std::string &field_function, const double target, const std::vector< unsigned int > &wires, const bool print=true)
 
void SetOptimisationParameters (const double dist=1., const double eps=1.e-4, const unsigned int nMaxIter=10)
 
bool ElectricFieldAtWire (const unsigned int iw, double &ex, double &ey)
 
void SetScanningGrid (const unsigned int nX, const unsigned int nY)
 
void SetScanningArea (const double xmin, const double xmax, const double ymin, const double ymax)
 
void SetScanningAreaLargest ()
 
void SetScanningAreaFirstOrder (const double scale=2.)
 
void EnableExtrapolation (const bool on=true)
 
void SetGravity (const double dx, const double dy, const double dz)
 Set the gravity orientation.
 
void GetGravity (double &dx, double &dy, double &dz) const
 Get the gravity orientation.
 
void EnableGravity (const bool on=true)
 Include gravity in the sag computation or not.
 
void EnableElectrostaticForce (const bool on=true)
 Include the electrostatic force in the sag computation or not.
 
bool ForcesOnWire (const unsigned int iw, std::vector< double > &xMap, std::vector< double > &yMap, std::vector< std::vector< double > > &fxMap, std::vector< std::vector< double > > &fyMap)
 
bool WireDisplacement (const unsigned int iw, const bool detailed, std::vector< double > &csag, std::vector< double > &xsag, std::vector< double > &ysag, double &stretch, const bool print=true)
 
void SetNumberOfShots (const unsigned int n)
 
void SetNumberOfSteps (const unsigned int n)
 Set the number of integration steps within each shot (must be >= 1).
 
Medium * GetMedium (const double x, const double y, const double z) override
 Get the medium at a given location (x, y, z).
 
void ElectricField (const double x, const double y, const double z, double &ex, double &ey, double &ez, Medium *&m, int &status) override
 
void ElectricField (const double x, const double y, const double z, double &ex, double &ey, double &ez, double &v, Medium *&m, int &status) override
 Calculate the drift field [V/cm] and potential [V] at (x, y, z).
 
bool GetVoltageRange (double &pmin, double &pmax) override
 Calculate the voltage range [V].
 
void WeightingField (const double x, const double y, const double z, double &wx, double &wy, double &wz, const std::string &label) override
 
double WeightingPotential (const double x, const double y, const double z, const std::string &label) override
 
bool GetBoundingBox (double &x0, double &y0, double &z0, double &x1, double &y1, double &z1) override
 Get the bounding box coordinates.
 
bool GetElementaryCell (double &x0, double &y0, double &z0, double &x1, double &y1, double &z1) override
 Get the coordinates of the elementary cell.
 
bool CrossedWire (const double x0, const double y0, const double z0, const double x1, const double y1, const double z1, double &xc, double &yc, double &zc, const bool centre, double &rc) override
 
bool InTrapRadius (const double q0, const double x0, const double y0, const double z0, double &xw, double &yx, double &rw) override
 
bool CrossedPlane (const double x0, const double y0, const double z0, const double x1, const double y1, const double z1, double &xc, double &yc, double &zc) override
 
double StepSizeHint () override
 
std::array< double, 3 > ElectricField (const double x, const double y, const double z)
 Calculate the drift field [V/cm] at (x, y, z).
 
- Public Member Functions inherited from Garfield::Component
 Component ()=delete
 Default constructor.
 
 Component (const std::string &name)
 Constructor.
 
virtual ~Component ()
 Destructor.
 
virtual void SetGeometry (Geometry *geo)
 Define the geometry.
 
virtual void Clear ()
 Reset.
 
std::array< double, 3 > ElectricField (const double x, const double y, const double z)
 Calculate the drift field [V/cm] at (x, y, z).
 
virtual double ElectricPotential (const double x, const double y, const double z)
 Calculate the (drift) electrostatic potential [V] at (x, y, z).
 
virtual void DelayedWeightingField (const double x, const double y, const double z, const double t, double &wx, double &wy, double &wz, const std::string &label)
 
virtual double DelayedWeightingPotential (const double x, const double y, const double z, const double t, const std::string &label)
 
virtual void MagneticField (const double x, const double y, const double z, double &bx, double &by, double &bz, int &status)
 
void SetMagneticField (const double bx, const double by, const double bz)
 Set a constant magnetic field.
 
virtual bool IsReady ()
 Ready for use?
 
double IntegrateFluxCircle (const double xc, const double yc, const double r, const unsigned int nI=50)
 
double IntegrateFluxSphere (const double xc, const double yc, const double zc, const double r, const unsigned int nI=20)
 
double IntegrateFluxParallelogram (const double x0, const double y0, const double z0, const double dx1, const double dy1, const double dz1, const double dx2, const double dy2, const double dz2, const unsigned int nU=20, const unsigned int nV=20)
 
double IntegrateWeightingFluxParallelogram (const std::string &label, const double x0, const double y0, const double z0, const double dx1, const double dy1, const double dz1, const double dx2, const double dy2, const double dz2, const unsigned int nU=20, const unsigned int nV=20)
 
double IntegrateFluxLine (const double x0, const double y0, const double z0, const double x1, const double y1, const double z1, const double xp, const double yp, const double zp, const unsigned int nI, const int isign=0)
 
void EnablePeriodicityX (const bool on=true)
 Enable simple periodicity in the $x$ direction.
 
void EnablePeriodicityY (const bool on=true)
 Enable simple periodicity in the $y$ direction.
 
void EnablePeriodicityZ (const bool on=true)
 Enable simple periodicity in the $z$ direction.
 
void IsPeriodic (bool &perx, bool &pery, bool &perz)
 Return periodicity flags.
 
void EnableMirrorPeriodicityX (const bool on=true)
 Enable mirror periodicity in the $x$ direction.
 
void EnableMirrorPeriodicityY (const bool on=true)
 Enable mirror periodicity in the $y$ direction.
 
void EnableMirrorPeriodicityZ (const bool on=true)
 Enable mirror periodicity in the $y$ direction.
 
void IsMirrorPeriodic (bool &perx, bool &pery, bool &perz)
 Return mirror periodicity flags.
 
void EnableAxialPeriodicityX (const bool on=true)
 Enable axial periodicity in the $x$ direction.
 
void EnableAxialPeriodicityY (const bool on=true)
 Enable axial periodicity in the $y$ direction.
 
void EnableAxialPeriodicityZ (const bool on=true)
 Enable axial periodicity in the $z$ direction.
 
void IsAxiallyPeriodic (bool &perx, bool &pery, bool &perz)
 Return axial periodicity flags.
 
void EnableRotationSymmetryX (const bool on=true)
 Enable rotation symmetry around the $x$ axis.
 
void EnableRotationSymmetryY (const bool on=true)
 Enable rotation symmetry around the $y$ axis.
 
void EnableRotationSymmetryZ (const bool on=true)
 Enable rotation symmetry around the $z$ axis.
 
void IsRotationSymmetric (bool &rotx, bool &roty, bool &rotz)
 Return rotation symmetry flags.
 
void EnableDebugging ()
 Switch on debugging messages.
 
void DisableDebugging ()
 Switch off debugging messages.
 
virtual bool HasMagneticField () const
 Does the component have a non-zero magnetic field?
 
virtual bool HasTownsendMap () const
 Does the component have maps of the Townsend coefficient?
 
virtual bool HasAttachmentMap () const
 Does the component have attachment maps?
 
virtual bool HasVelocityMap () const
 Does the component have velocity maps?
 
virtual bool ElectronAttachment (const double, const double, const double, double &eta)
 Get the electron attachment coefficient.
 
virtual bool HoleAttachment (const double, const double, const double, double &eta)
 Get the hole attachment coefficient.
 
virtual bool ElectronTownsend (const double, const double, const double, double &alpha)
 Get the electron Townsend coefficient.
 
virtual bool HoleTownsend (const double, const double, const double, double &alpha)
 Get the hole Townsend coefficient.
 
virtual bool ElectronVelocity (const double, const double, const double, double &vx, double &vy, double &vz)
 Get the electron drift velocity.
 
virtual bool HoleVelocity (const double, const double, const double, double &vx, double &vy, double &vz)
 Get the hole drift velocity.
 
virtual bool GetElectronLifetime (const double, const double, const double, double &etau)
 
virtual bool GetHoleLifetime (const double, const double, const double, double &htau)
 

Additional Inherited Members

- Protected Attributes inherited from Garfield::Component
std::string m_className = "Component"
 Class name.
 
Geometry * m_geometry = nullptr
 Pointer to the geometry.
 
std::array< double, 3 > m_b0 = {{0., 0., 0.}}
 Constant magnetic field.
 
bool m_ready = false
 Ready for use?
 
bool m_debug = false
 Switch on/off debugging messages.
 
std::array< bool, 3 > m_periodic = {{false, false, false}}
 Simple periodicity in x, y, z.
 
std::array< bool, 3 > m_mirrorPeriodic = {{false, false, false}}
 Mirror periodicity in x, y, z.
 
std::array< bool, 3 > m_axiallyPeriodic = {{false, false, false}}
 Axial periodicity in x, y, z.
 
std::array< bool, 3 > m_rotationSymmetric = {{false, false, false}}
 Rotation symmetry around x-axis, y-axis, z-axis.
 

Detailed Description

Semi-analytic calculation of two-dimensional configurations consisting of wires, planes, and tubes.

Definition at line 18 of file ComponentAnalyticField.hh.

Member Enumeration Documentation

◆ Cell

enum Garfield::ComponentAnalyticField::Cell
Enumerator
A00 
B1X 
B1Y 
B2X 
B2Y 
C10 
C2X 
C2Y 
C30 
D10 
D20 
D30 
D40 
Unknown 

Definition at line 392 of file ComponentAnalyticField.hh.

392 {
393 A00,
394 B1X,
395 B1Y,
396 B2X,
397 B2Y,
398 C10,
399 C2X,
400 C2Y,
401 C30,
402 D10,
403 D20,
404 D30,
405 D40,
406 Unknown
407 };

Constructor & Destructor Documentation

◆ ComponentAnalyticField()

Garfield::ComponentAnalyticField::ComponentAnalyticField ( )

Constructor.

Definition at line 279 of file ComponentAnalyticField.cc.

279 : Component("AnalyticField") {
280 CellInit();
281}
Garfield::Component::Component
Component()=delete
Default constructor.

◆ ~ComponentAnalyticField()

Garfield::ComponentAnalyticField::~ComponentAnalyticField ( )
inline

Destructor.

Definition at line 23 of file ComponentAnalyticField.hh.

23{}

Member Function Documentation

◆ AddCharge()

void Garfield::ComponentAnalyticField::AddCharge ( const double x,
const double y,
const double z,
const double q )

Add a point charge.

Definition at line 1764 of file ComponentAnalyticField.cc.

1765 {
1766 // Convert from fC to internal units (division by 4 pi epsilon0).
1767 Charge3d charge;
1768 charge.x = x;
1769 charge.y = y;
1770 charge.z = z;
1771 charge.e = q / FourPiEpsilon0;
1772 m_ch3d.push_back(std::move(charge));
1773}

◆ AddPixelOnPlanePhi()

void Garfield::ComponentAnalyticField::AddPixelOnPlanePhi ( const double phi,
const double rmin,
const double rmax,
const double zmin,
const double zmax,
const std::string & label,
const double gap = -1. )

Add a pixel on an existing plane at constant phi.

Definition at line 1547 of file ComponentAnalyticField.cc.

1550 {
1551 if (!m_polar || (!m_ynplan[2] && !m_ynplan[3])) {
1552 std::cerr << m_className << "::AddPixelOnPlanePhi:\n"
1553 << " There are no planes at constant phi.\n";
1554 return;
1555 }
1556
1557 if (fabs(rmax - rmin) < Small || fabs(zmax - zmin) < Small) {
1558 std::cerr << m_className << "::AddPixelOnPlanePhi:\n"
1559 << " Pixel width must be greater than zero.\n";
1560 return;
1561 }
1562 if (rmin < Small || rmax < Small) {
1563 std::cerr << m_className << "::AddPixelOnPlanePhi:\n"
1564 << " Radius must be greater than zero.\n";
1565 return;
1566 }
1567
1568 if (label.empty()) {
1569 std::cerr << m_className << "::AddPixelOnPlanePhi: "
1570 << "Label must not be an empty string.\n";
1571 return;
1572 }
1573
1574 Pixel pixel;
1575 pixel.type = label;
1576 pixel.ind = -1;
1577 const double smin = log(rmin);
1578 const double smax = log(rmax);
1579 pixel.smin = std::min(smin, smax);
1580 pixel.smax = std::max(smin, smax);
1581 pixel.zmin = std::min(zmin, zmax);
1582 pixel.zmax = std::max(zmin, zmax);
1583 pixel.gap = gap > Small ? DegreeToRad * gap : -1.;
1584
1585 int iplane = 2;
1586 if (m_ynplan[3]) {
1587 const double d0 = fabs(m_coplan[2] - phi * DegreeToRad);
1588 const double d1 = fabs(m_coplan[3] - phi * DegreeToRad);
1589 if (d1 < d0) iplane = 3;
1590 }
1591
1592 m_planes[iplane].pixels.push_back(std::move(pixel));
1593 // Add the identifier to the list of readout groups.
1594 AddReadout(label, true);
1595}
Garfield::ComponentAnalyticField::AddReadout
void AddReadout(const std::string &label, const bool silent=false)
Definition ComponentAnalyticField.cc:3772
Garfield::Component::m_className
std::string m_className
Class name.
Definition Component.hh:359
Heed::fabs
DoubleAc fabs(const DoubleAc &f)
Definition DoubleAc.h:615

◆ AddPixelOnPlaneR()

void Garfield::ComponentAnalyticField::AddPixelOnPlaneR ( const double r,
const double phimin,
const double phimax,
const double zmin,
const double zmax,
const std::string & label,
const double gap = -1. )

Add a pixel on an existing plane at constant radius.

Definition at line 1501 of file ComponentAnalyticField.cc.

1504 {
1505 if (!m_polar || (!m_ynplan[0] && !m_ynplan[1])) {
1506 std::cerr << m_className << "::AddPixelOnPlaneR:\n"
1507 << " There are no planes at constant r.\n";
1508 return;
1509 }
1510
1511 if (fabs(phimax - phimin) < Small || fabs(zmax - zmin) < Small) {
1512 std::cerr << m_className << "::AddPixelOnPlaneR:\n"
1513 << " Pixel width must be greater than zero.\n";
1514 return;
1515 }
1516
1517 if (label.empty()) {
1518 std::cerr << m_className << "::AddPixelOnPlaneR: "
1519 << "Label must not be an empty string.\n";
1520 return;
1521 }
1522
1523 Pixel pixel;
1524 pixel.type = label;
1525 pixel.ind = -1;
1526 const double smin = phimin * DegreeToRad;
1527 const double smax = phimax * DegreeToRad;
1528 pixel.smin = std::min(smin, smax);
1529 pixel.smax = std::max(smin, smax);
1530 pixel.zmin = std::min(zmin, zmax);
1531 pixel.zmax = std::max(zmin, zmax);
1532 pixel.gap = gap > Small ? gap : -1.;
1533
1534 int iplane = 0;
1535 if (m_ynplan[1]) {
1536 const double rho = r > 0. ? log(r) : -25.;
1537 const double d0 = fabs(m_coplan[0] - rho);
1538 const double d1 = fabs(m_coplan[1] - rho);
1539 if (d1 < d0) iplane = 1;
1540 }
1541
1542 m_planes[iplane].pixels.push_back(std::move(pixel));
1543 // Add the identifier to the list of readout groups.
1544 AddReadout(label, true);
1545}

◆ AddPixelOnPlaneX()

void Garfield::ComponentAnalyticField::AddPixelOnPlaneX ( const double x,
const double ymin,
const double ymax,
const double zmin,
const double zmax,
const std::string & label,
const double gap = -1.,
const double rot = 0. )

Add a pixel on an existing plane at constant x.

Parameters
xcoordinate of the plane.
yminlower limit of the pixel cell in y,
ymaxupper limit of the pixel cell in y.
zminlower limit of the pixel cell in z.
zmaxupper limit of the pixel cell in z.
labelweighting field identifier.
gapdistance to the opposite plane (optional).
rotrotation angle (rad) of the pixel (optional).

Definition at line 1407 of file ComponentAnalyticField.cc.

1410 {
1411 if (m_polar || (!m_ynplan[0] && !m_ynplan[1])) {
1412 std::cerr << m_className << "::AddPixelOnPlaneX:\n"
1413 << " There are no planes at constant x.\n";
1414 return;
1415 }
1416
1417 if (fabs(ymax - ymin) < Small || fabs(zmax - zmin) < Small) {
1418 std::cerr << m_className << "::AddPixelOnPlaneX:\n"
1419 << " Pixel width must be greater than zero.\n";
1420 return;
1421 }
1422
1423 if (label.empty()) {
1424 std::cerr << m_className << "::AddPixelOnPlaneX: "
1425 << "Label must not be an empty string.\n";
1426 return;
1427 }
1428
1429 Pixel pixel;
1430 pixel.type = label;
1431 pixel.ind = -1;
1432 pixel.smin = std::min(ymin, ymax);
1433 pixel.smax = std::max(ymin, ymax);
1434 pixel.zmin = std::min(zmin, zmax);
1435 pixel.zmax = std::max(zmin, zmax);
1436 pixel.gap = gap > Small ? gap : -1.;
1437 if (fabs(rot) > 1.e-9) {
1438 pixel.cphi = cos(rot);
1439 pixel.sphi = sin(rot);
1440 }
1441
1442 int iplane = 0;
1443 if (m_ynplan[1]) {
1444 const double d0 = fabs(m_coplan[0] - x);
1445 const double d1 = fabs(m_coplan[1] - x);
1446 if (d1 < d0) iplane = 1;
1447 }
1448
1449 m_planes[iplane].pixels.push_back(std::move(pixel));
1450 // Add the identifier to the list of readout groups.
1451 AddReadout(label, true);
1452}
Heed::cos
DoubleAc cos(const DoubleAc &f)
Definition DoubleAc.cpp:432
Heed::sin
DoubleAc sin(const DoubleAc &f)
Definition DoubleAc.cpp:384

◆ AddPixelOnPlaneY()

void Garfield::ComponentAnalyticField::AddPixelOnPlaneY ( const double y,
const double xmin,
const double xmax,
const double zmin,
const double zmax,
const std::string & label,
const double gap = -1.,
const double rot = 0. )

Add a pixel on an existing plane at constant y.

Definition at line 1454 of file ComponentAnalyticField.cc.

1457 {
1458 if (m_polar || (!m_ynplan[2] && !m_ynplan[3])) {
1459 std::cerr << m_className << "::AddPixelOnPlaneY:\n"
1460 << " There are no planes at constant y.\n";
1461 return;
1462 }
1463
1464 if (fabs(xmax - xmin) < Small || fabs(zmax - zmin) < Small) {
1465 std::cerr << m_className << "::AddPixelOnPlaneY:\n"
1466 << " Pixel width must be greater than zero.\n";
1467 return;
1468 }
1469
1470 if (label.empty()) {
1471 std::cerr << m_className << "::AddPixelOnPlaneY: "
1472 << "Label must not be an empty string.\n";
1473 return;
1474 }
1475
1476 Pixel pixel;
1477 pixel.type = label;
1478 pixel.ind = -1;
1479 pixel.smin = std::min(xmin, xmax);
1480 pixel.smax = std::max(xmin, xmax);
1481 pixel.zmin = std::min(zmin, zmax);
1482 pixel.zmax = std::max(zmin, zmax);
1483 pixel.gap = gap > Small ? gap : -1.;
1484 if (fabs(rot) > 1.e-9) {
1485 pixel.cphi = cos(rot);
1486 pixel.sphi = sin(rot);
1487 }
1488
1489 int iplane = 2;
1490 if (m_ynplan[3]) {
1491 const double d0 = fabs(m_coplan[2] - y);
1492 const double d1 = fabs(m_coplan[3] - y);
1493 if (d1 < d0) iplane = 3;
1494 }
1495
1496 m_planes[iplane].pixels.push_back(std::move(pixel));
1497 // Add the identifier to the list of readout groups.
1498 AddReadout(label, true);
1499}

◆ AddPlanePhi()

void Garfield::ComponentAnalyticField::AddPlanePhi ( const double phi,
const double voltage,
const std::string & label = "" )

Add a plane at constant phi.

Definition at line 1135 of file ComponentAnalyticField.cc.

1136 {
1137 if (!m_polar) {
1138 std::cerr << m_className << "::AddPlanePhi:\n"
1139 << " Not compatible with Cartesian coordinates; ignored.\n";
1140 return;
1141 }
1142 if (m_ynplan[2] && m_ynplan[3]) {
1143 std::cerr << m_className << "::AddPlanePhi:\n"
1144 << " Cannot have more than two planes at constant phi.\n";
1145 return;
1146 }
1147
1148 if (m_ynplan[2]) {
1149 m_ynplan[3] = true;
1150 m_coplan[3] = phi * DegreeToRad;
1151 m_vtplan[3] = v;
1152 m_planes[3].type = label;
1153 m_planes[3].ind = -1;
1154 } else {
1155 m_ynplan[2] = true;
1156 m_coplan[2] = phi * DegreeToRad;
1157 m_vtplan[2] = v;
1158 m_planes[2].type = label;
1159 m_planes[2].ind = -1;
1160 // Switch off default periodicity.
1161 if (m_pery && std::abs(m_sy - TwoPi) < 1.e-4) {
1162 m_pery = false;
1163 }
1164 }
1165
1166 // Add the identifier to the list of readout groups.
1167 if (!label.empty()) AddReadout(label, true);
1168 // Force recalculation of the capacitance and signal matrices.
1169 m_cellset = false;
1170 m_sigset = false;
1171}

◆ AddPlaneR()

void Garfield::ComponentAnalyticField::AddPlaneR ( const double r,
const double voltage,
const std::string & label = "" )

Add a plane at constant radius.

Definition at line 1095 of file ComponentAnalyticField.cc.

1096 {
1097 if (!m_polar) {
1098 std::cerr << m_className << "::AddPlaneR:\n"
1099 << " Not compatible with Cartesian coordinates; ignored.\n";
1100 return;
1101 }
1102 if (r <= 0.) {
1103 std::cerr << m_className << "::AddPlaneR:\n"
1104 << " Radius must be larger than zero; plane ignored.\n";
1105 return;
1106 }
1107
1108 if (m_ynplan[0] && m_ynplan[1]) {
1109 std::cerr << m_className << "::AddPlaneR:\n"
1110 << " Cannot have more than two circular planes.\n";
1111 return;
1112 }
1113
1114 if (m_ynplan[0]) {
1115 m_ynplan[1] = true;
1116 m_coplan[1] = log(r);
1117 m_vtplan[1] = v;
1118 m_planes[1].type = label;
1119 m_planes[1].ind = -1;
1120 } else {
1121 m_ynplan[0] = true;
1122 m_coplan[0] = log(r);
1123 m_vtplan[0] = v;
1124 m_planes[0].type = label;
1125 m_planes[0].ind = -1;
1126 }
1127
1128 // Add the identifier to the list of readout groups.
1129 if (!label.empty()) AddReadout(label, true);
1130 // Force recalculation of the capacitance and signal matrices.
1131 m_cellset = false;
1132 m_sigset = false;
1133}

◆ AddPlaneX()

void Garfield::ComponentAnalyticField::AddPlaneX ( const double x,
const double voltage,
const std::string & label = "" )

Add a plane at constant x.

Definition at line 1027 of file ComponentAnalyticField.cc.

1028 {
1029 if (m_polar) {
1030 std::cerr << m_className << "::AddPlaneX:\n"
1031 << " Not compatible with polar coordinates; ignored.\n";
1032 return;
1033 }
1034 if (m_ynplan[0] && m_ynplan[1]) {
1035 std::cerr << m_className << "::AddPlaneX:\n"
1036 << " Cannot have more than two planes at constant x.\n";
1037 return;
1038 }
1039
1040 if (m_ynplan[0]) {
1041 m_ynplan[1] = true;
1042 m_coplan[1] = x;
1043 m_vtplan[1] = v;
1044 m_planes[1].type = label;
1045 m_planes[1].ind = -1;
1046 } else {
1047 m_ynplan[0] = true;
1048 m_coplan[0] = x;
1049 m_vtplan[0] = v;
1050 m_planes[0].type = label;
1051 m_planes[0].ind = -1;
1052 }
1053
1054 // Add the identifier to the list of readout groups.
1055 if (!label.empty()) AddReadout(label, true);
1056 // Force recalculation of the capacitance and signal matrices.
1057 m_cellset = false;
1058 m_sigset = false;
1059}

◆ AddPlaneY()

void Garfield::ComponentAnalyticField::AddPlaneY ( const double y,
const double voltage,
const std::string & label = "" )

Add a plane at constant y.

Definition at line 1061 of file ComponentAnalyticField.cc.

1062 {
1063 if (m_polar) {
1064 std::cerr << m_className << "::AddPlaneY:\n"
1065 << " Not compatible with polar coordinates; ignored.\n";
1066 return;
1067 }
1068 if (m_ynplan[2] && m_ynplan[3]) {
1069 std::cerr << m_className << "::AddPlaneY:\n"
1070 << " Cannot have more than two planes at constant y.\n";
1071 return;
1072 }
1073
1074 if (m_ynplan[2]) {
1075 m_ynplan[3] = true;
1076 m_coplan[3] = y;
1077 m_vtplan[3] = v;
1078 m_planes[3].type = label;
1079 m_planes[3].ind = -1;
1080 } else {
1081 m_ynplan[2] = true;
1082 m_coplan[2] = y;
1083 m_vtplan[2] = v;
1084 m_planes[2].type = label;
1085 m_planes[2].ind = -1;
1086 }
1087
1088 // Add the identifier to the list of readout groups.
1089 if (!label.empty()) AddReadout(label, true);
1090 // Force recalculation of the capacitance and signal matrices.
1091 m_cellset = false;
1092 m_sigset = false;
1093}

Referenced by main().

◆ AddReadout()

void Garfield::ComponentAnalyticField::AddReadout ( const std::string & label,
const bool silent = false )

Request calculation of weighting field and potential for a given group of wires or planes.

Definition at line 3772 of file ComponentAnalyticField.cc.

3773 {
3774 // Check if this readout group already exists.
3775 if (std::find(m_readout.begin(), m_readout.end(), label) != m_readout.end()) {
3776 if (!silent) {
3777 std::cout << m_className << "::AddReadout: Readout group "
3778 << label << " already exists.\n";
3779 }
3780 } else {
3781 m_readout.push_back(label);
3782 }
3783
3784 unsigned int nWiresFound = 0;
3785 for (const auto& wire : m_w) {
3786 if (wire.type == label) ++nWiresFound;
3787 }
3788
3789 unsigned int nPlanesFound = 0;
3790 unsigned int nStripsFound = 0;
3791 unsigned int nPixelsFound = 0;
3792 for (const auto& plane : m_planes) {
3793 if (plane.type == label) ++nPlanesFound;
3794 for (const auto& strip : plane.strips1) {
3795 if (strip.type == label) ++nStripsFound;
3796 }
3797 for (const auto& strip : plane.strips2) {
3798 if (strip.type == label) ++nStripsFound;
3799 }
3800 for (const auto& pixel : plane.pixels) {
3801 if (pixel.type == label) ++nPixelsFound;
3802 }
3803 }
3804
3805 if (nWiresFound == 0 && nPlanesFound == 0 && nStripsFound == 0 &&
3806 nPixelsFound == 0) {
3807 std::cerr << m_className << "::AddReadout:\n"
3808 << " At present there are no wires, planes or strips\n"
3809 << " associated to readout group " << label << ".\n";
3810 } else if (!silent) {
3811 std::cout << m_className << "::AddReadout:\n"
3812 << " Readout group " << label << " comprises:\n";
3813 if (nWiresFound > 1) {
3814 std::cout << " " << nWiresFound << " wires\n";
3815 } else if (nWiresFound == 1) {
3816 std::cout << " 1 wire\n";
3817 }
3818 if (nPlanesFound > 1) {
3819 std::cout << " " << nPlanesFound << " planes\n";
3820 } else if (nPlanesFound == 1) {
3821 std::cout << " 1 plane\n";
3822 }
3823 if (nStripsFound > 1) {
3824 std::cout << " " << nStripsFound << " strips\n";
3825 } else if (nStripsFound == 1) {
3826 std::cout << " 1 strip\n";
3827 }
3828 if (nPixelsFound > 1) {
3829 std::cout << " " << nPixelsFound << " pixels\n";
3830 } else if (nPixelsFound == 1) {
3831 std::cout << " 1 pixel\n";
3832 }
3833 }
3834 m_sigset = false;
3835}

Referenced by AddPixelOnPlanePhi(), AddPixelOnPlaneR(), AddPixelOnPlaneX(), AddPixelOnPlaneY(), AddPlanePhi(), AddPlaneR(), AddPlaneX(), AddPlaneY(), AddStripOnPlanePhi(), AddStripOnPlaneR(), AddStripOnPlaneX(), AddStripOnPlaneY(), AddTube(), and AddWire().

◆ AddStripOnPlanePhi()

void Garfield::ComponentAnalyticField::AddStripOnPlanePhi ( const char direction,
const double phi,
const double smin,
const double smax,
const std::string & label,
const double gap = -1. )

Add a strip in the r or z direction on an existing plane at constant phi.

Definition at line 1340 of file ComponentAnalyticField.cc.

1345 {
1346 if (!m_polar || (!m_ynplan[2] && !m_ynplan[3])) {
1347 std::cerr << m_className << "::AddStripOnPlanePhi:\n"
1348 << " There are no planes at constant phi.\n";
1349 return;
1350 }
1351
1352 if (dir != 'r' && dir != 'R' && dir != 'z' && dir != 'Z') {
1353 std::cerr << m_className << "::AddStripOnPlanePhi:\n"
1354 << " Invalid direction (" << dir << ").\n"
1355 << " Only strips in r or z direction are possible.\n";
1356 return;
1357 }
1358
1359 if (fabs(smax - smin) < Small) {
1360 std::cerr << m_className << "::AddStripOnPlanePhi:\n"
1361 << " Strip width must be greater than zero.\n";
1362 return;
1363 }
1364
1365 if (label.empty()) {
1366 std::cerr << m_className << "::AddStripOnPlanePhi: "
1367 << "Label must not be an empty string.\n";
1368 return;
1369 }
1370
1371 Strip newStrip;
1372 newStrip.type = label;
1373 newStrip.ind = -1;
1374 if (dir== 'z' || dir == 'Z') {
1375 if (smin < Small || smax < Small) {
1376 std::cerr << m_className << "::AddStripOnPlanePhi:\n"
1377 << " Radius must be greater than zero.\n";
1378 return;
1379 }
1380 const double rhomin = log(smin);
1381 const double rhomax = log(smax);
1382 newStrip.smin = std::min(rhomin, rhomax);
1383 newStrip.smax = std::max(rhomin, rhomax);
1384 } else {
1385 newStrip.smin = std::min(smin, smax);
1386 newStrip.smax = std::max(smin, smax);
1387 }
1388 newStrip.gap = gap > Small ? DegreeToRad * gap : -1.;
1389
1390 int iplane = 2;
1391 if (m_ynplan[3]) {
1392 const double d2 = fabs(m_coplan[2] - phi * DegreeToRad);
1393 const double d3 = fabs(m_coplan[3] - phi * DegreeToRad);
1394 if (d3 < d2) iplane = 3;
1395 }
1396
1397 if (dir == 'r' || dir == 'R') {
1398 m_planes[iplane].strips1.push_back(std::move(newStrip));
1399 } else {
1400 m_planes[iplane].strips2.push_back(std::move(newStrip));
1401 }
1402 // Add the identifier to the list of readout groups.
1403 AddReadout(label, true);
1404}

◆ AddStripOnPlaneR()

void Garfield::ComponentAnalyticField::AddStripOnPlaneR ( const char direction,
const double r,
const double smin,
const double smax,
const std::string & label,
const double gap = -1. )

Add a strip in the phi or z direction on an existing plane at constant radius.

Definition at line 1279 of file ComponentAnalyticField.cc.

1283 {
1284 if (!m_polar || (!m_ynplan[0] && !m_ynplan[1])) {
1285 std::cerr << m_className << "::AddStripOnPlaneR:\n"
1286 << " There are no planes at constant r.\n";
1287 return;
1288 }
1289
1290 if (dir != 'p' && dir != 'P' && dir != 'z' && dir != 'Z') {
1291 std::cerr << m_className << "::AddStripOnPlaneR:\n"
1292 << " Invalid direction (" << dir << ").\n"
1293 << " Only strips in p(hi) or z direction are possible.\n";
1294 return;
1295 }
1296
1297 if (fabs(smax - smin) < Small) {
1298 std::cerr << m_className << "::AddStripOnPlaneR:\n"
1299 << " Strip width must be greater than zero.\n";
1300 return;
1301 }
1302
1303 if (label.empty()) {
1304 std::cerr << m_className << "::AddStripOnPlaneR: "
1305 << "Label must not be an empty string.\n";
1306 return;
1307 }
1308
1309 Strip newStrip;
1310 newStrip.type = label;
1311 newStrip.ind = -1;
1312 if (dir == 'z' || dir == 'Z') {
1313 const double phimin = smin * DegreeToRad;
1314 const double phimax = smax * DegreeToRad;
1315 newStrip.smin = std::min(phimin, phimax);
1316 newStrip.smax = std::max(phimin, phimax);
1317 } else {
1318 newStrip.smin = std::min(smin, smax);
1319 newStrip.smax = std::max(smin, smax);
1320 }
1321 newStrip.gap = gap > Small ? gap : -1.;
1322
1323 int iplane = 0;
1324 if (m_ynplan[1]) {
1325 const double rho = r > 0. ? log(r) : -25.;
1326 const double d0 = fabs(m_coplan[0] - rho);
1327 const double d1 = fabs(m_coplan[1] - rho);
1328 if (d1 < d0) iplane = 1;
1329 }
1330
1331 if (dir == 'p' || dir == 'P') {
1332 m_planes[iplane].strips1.push_back(std::move(newStrip));
1333 } else {
1334 m_planes[iplane].strips2.push_back(std::move(newStrip));
1335 }
1336 // Add the identifier to the list of readout groups.
1337 AddReadout(label, true);
1338}

◆ AddStripOnPlaneX()

void Garfield::ComponentAnalyticField::AddStripOnPlaneX ( const char direction,
const double x,
const double smin,
const double smax,
const std::string & label,
const double gap = -1. )

Add a strip in the y or z direction on an existing plane at constant x.

Parameters
direction'y' or 'z'.
xcoordinate of the plane.
sminlower limit of the strip in y or z.
smaxupper limit of the strip in y or z.
labelweighting field identifier.
gapdistance to the opposite plane (optional).

Definition at line 1173 of file ComponentAnalyticField.cc.

1177 {
1178 if (m_polar || (!m_ynplan[0] && !m_ynplan[1])) {
1179 std::cerr << m_className << "::AddStripOnPlaneX:\n"
1180 << " There are no planes at constant x.\n";
1181 return;
1182 }
1183
1184 if (dir != 'y' && dir != 'Y' && dir != 'z' && dir != 'Z') {
1185 std::cerr << m_className << "::AddStripOnPlaneX:\n"
1186 << " Invalid direction (" << dir << ").\n"
1187 << " Only strips in y or z direction are possible.\n";
1188 return;
1189 }
1190
1191 if (fabs(smax - smin) < Small) {
1192 std::cerr << m_className << "::AddStripOnPlaneX:\n"
1193 << " Strip width must be greater than zero.\n";
1194 return;
1195 }
1196
1197 if (label.empty()) {
1198 std::cerr << m_className << "::AddStripOnPlaneX: "
1199 << "Label must not be an empty string.\n";
1200 return;
1201 }
1202
1203 Strip newStrip;
1204 newStrip.type = label;
1205 newStrip.ind = -1;
1206 newStrip.smin = std::min(smin, smax);
1207 newStrip.smax = std::max(smin, smax);
1208 newStrip.gap = gap > Small ? gap : -1.;
1209
1210 int iplane = 0;
1211 if (m_ynplan[1]) {
1212 const double d0 = fabs(m_coplan[0] - x);
1213 const double d1 = fabs(m_coplan[1] - x);
1214 if (d1 < d0) iplane = 1;
1215 }
1216
1217 if (dir == 'y' || dir == 'Y') {
1218 m_planes[iplane].strips1.push_back(std::move(newStrip));
1219 } else {
1220 m_planes[iplane].strips2.push_back(std::move(newStrip));
1221 }
1222 // Add the identifier to the list of readout groups.
1223 AddReadout(label, true);
1224}

◆ AddStripOnPlaneY()

void Garfield::ComponentAnalyticField::AddStripOnPlaneY ( const char direction,
const double y,
const double smin,
const double smax,
const std::string & label,
const double gap = -1. )

Add a strip in the x or z direction on an existing plane at constant y.

Definition at line 1226 of file ComponentAnalyticField.cc.

1230 {
1231 if (m_polar || (!m_ynplan[2] && !m_ynplan[3])) {
1232 std::cerr << m_className << "::AddStripOnPlaneY:\n"
1233 << " There are no planes at constant y.\n";
1234 return;
1235 }
1236
1237 if (dir != 'x' && dir != 'X' && dir != 'z' && dir != 'Z') {
1238 std::cerr << m_className << "::AddStripOnPlaneY:\n"
1239 << " Invalid direction (" << dir << ").\n"
1240 << " Only strips in x or z direction are possible.\n";
1241 return;
1242 }
1243
1244 if (fabs(smax - smin) < Small) {
1245 std::cerr << m_className << "::AddStripOnPlaneY:\n"
1246 << " Strip width must be greater than zero.\n";
1247 return;
1248 }
1249
1250 if (label.empty()) {
1251 std::cerr << m_className << "::AddStripOnPlaneY: "
1252 << "Label must not be an empty string.\n";
1253 return;
1254 }
1255
1256 Strip newStrip;
1257 newStrip.type = label;
1258 newStrip.ind = -1;
1259 newStrip.smin = std::min(smin, smax);
1260 newStrip.smax = std::max(smin, smax);
1261 newStrip.gap = gap > Small ? gap : -1.;
1262
1263 int iplane = 2;
1264 if (m_ynplan[3]) {
1265 const double d2 = fabs(m_coplan[2] - y);
1266 const double d3 = fabs(m_coplan[3] - y);
1267 if (d3 < d2) iplane = 3;
1268 }
1269
1270 if (dir == 'x' || dir == 'X') {
1271 m_planes[iplane].strips1.push_back(std::move(newStrip));
1272 } else {
1273 m_planes[iplane].strips2.push_back(std::move(newStrip));
1274 }
1275 // Add the identifier to the list of readout groups.
1276 AddReadout(label, true);
1277}

Referenced by main().

◆ AddTube()

void Garfield::ComponentAnalyticField::AddTube ( const double radius,
const double voltage,
const int nEdges,
const std::string & label = "" )

Add a tube.

Definition at line 985 of file ComponentAnalyticField.cc.

987 {
988 // Check if the provided parameters make sense.
989 if (radius <= 0.0) {
990 std::cerr << m_className << "::AddTube: Unphysical tube dimension.\n";
991 return;
992 }
993
994 if (nEdges < 3 && nEdges != 0) {
995 std::cerr << m_className << "::AddTube: Unphysical number of tube edges ("
996 << nEdges << ")\n";
997 return;
998 }
999
1000 // If there is already a tube defined, print a warning message.
1001 if (m_tube) {
1002 std::cout << m_className << "::AddTube:\n"
1003 << " Warning: Existing tube settings will be overwritten.\n";
1004 }
1005
1006 // Set the coordinate system.
1007 m_tube = true;
1008 m_polar = false;
1009
1010 // Set the tube parameters.
1011 m_cotube = radius;
1012 m_cotube2 = radius * radius;
1013 m_vttube = voltage;
1014
1015 m_ntube = nEdges;
1016
1017 m_planes[4].type = label;
1018 m_planes[4].ind = -1;
1019
1020 // Add the identifier to the list of readout groups.
1021 if (!label.empty()) AddReadout(label, true);
1022 // Force recalculation of the capacitance and signal matrices.
1023 m_cellset = false;
1024 m_sigset = false;
1025}

◆ AddWire()

void Garfield::ComponentAnalyticField::AddWire ( const double x,
const double y,
const double diameter,
const double voltage,
const std::string & label = "",
const double length = 100.,
const double tension = 50.,
const double rho = 19.3,
const int ntrap = 5 )

Add a wire at (x, y) .

Definition at line 917 of file ComponentAnalyticField.cc.

922 {
923 // Check if the provided parameters make sense.
924 if (diameter <= 0.) {
925 std::cerr << m_className << "::AddWire: Unphysical wire diameter.\n";
926 return;
927 }
928
929 if (tension <= 0.) {
930 std::cerr << m_className << "::AddWire: Unphysical wire tension.\n";
931 return;
932 }
933
934 if (rho <= 0.) {
935 std::cerr << m_className << "::AddWire: Unphysical wire density.\n";
936 return;
937 }
938
939 if (length <= 0.) {
940 std::cerr << m_className << "::AddWire: Unphysical wire length.\n";
941 return;
942 }
943
944 if (ntrap <= 0) {
945 std::cerr << m_className << "::AddWire: Nbr. of trap radii must be > 0.\n";
946 return;
947 }
948
949 if (m_polar && x <= diameter) {
950 std::cerr << m_className << "::AddWire: Wire is too close to the origin.\n";
951 return;
952 }
953
954 // Create a new wire.
955 Wire wire;
956 if (m_polar) {
957 double r = 0., phi = 0.;
958 Polar2Internal(x, y, r, phi);
959 wire.x = r;
960 wire.y = phi;
961 wire.r = 0.5 * diameter / x;
962 } else {
963 wire.x = x;
964 wire.y = y;
965 wire.r = 0.5 * diameter;
966 }
967 wire.v = voltage;
968 wire.u = length;
969 wire.type = label;
970 wire.e = 0.;
971 wire.ind = -1;
972 wire.nTrap = ntrap;
973 wire.tension = tension;
974 wire.density = rho;
975 // Add the wire to the list.
976 m_w.push_back(std::move(wire));
977 ++m_nWires;
978 // Add the identifier to the list of readout groups.
979 if (!label.empty()) AddReadout(label, true);
980 // Force recalculation of the capacitance and signal matrices.
981 m_cellset = false;
982 m_sigset = false;
983}

◆ ClearCharges()

void Garfield::ComponentAnalyticField::ClearCharges ( )

Remove all point charges.

Definition at line 1775 of file ComponentAnalyticField.cc.

1775 {
1776 m_ch3d.clear();
1777 m_nTermBessel = 10;
1778 m_nTermPoly = 100;
1779}

◆ CrossedPlane()

bool Garfield::ComponentAnalyticField::CrossedPlane ( const double x0,
const double y0,
const double z0,
const double x1,
const double y1,
const double z1,
double & xc,
double & yc,
double & zc )
overridevirtual

Determine whether the line between two points crosses a plane.

Reimplemented from Garfield::Component.

Definition at line 878 of file ComponentAnalyticField.cc.

881 {
882
883 double x0 = xx0;
884 double y0 = yy0;
885 double x1 = xx1;
886 double y1 = yy1;
887 if (m_polar) {
888 Cartesian2Internal(xx0, yy0, x0, y0);
889 Cartesian2Internal(xx1, yy1, x1, y1);
890 }
891
892 double smin = -1.;
893 for (size_t i = 0; i < 2; ++i) {
894 if (!m_ynplan[i]) continue;
895 if (SameSide(x0, x1, m_coplan[i])) continue;
896 const double s = (m_coplan[i] - x0) / (x1 - x0);
897 if (smin < 0. || s < smin) {
898 smin = s;
899 }
900 }
901 for (size_t i = 2; i < 4; ++i) {
902 if (!m_ynplan[i]) continue;
903 if (SameSide(y0, y1, m_coplan[i])) continue;
904 const double s = (m_coplan[i] - y0) / (y1 - y0);
905 if (smin < 0. || s < smin) {
906 smin = s;
907 }
908 }
909 if (smin < 0.) return false;
910 xc = x0 + smin * (x1 - x0);
911 yc = y0 + smin * (y1 - y0);
912 zc = z0 + smin * (z1 - z0);
913 if (m_polar) Internal2Cartesian(xc, yc, xc, yc);
914 return true;
915}

◆ CrossedWire()

bool Garfield::ComponentAnalyticField::CrossedWire ( const double x0,
const double y0,
const double z0,
const double x1,
const double y1,
const double z1,
double & xc,
double & yc,
double & zc,
const bool centre,
double & rc )
overridevirtual

Determine whether the line between two points crosses a wire.

Parameters
x0,y0,z0first point [cm].
x1,y1,z1second point [cm]
xc,yc,zcpoint [cm] where the line crosses the wire or the coordinates of the wire centre.
centreflag whether to return the coordinates of the line-wire crossing point or of the wire centre.
rcradius [cm] of the wire.

Reimplemented from Garfield::Component.

Definition at line 698 of file ComponentAnalyticField.cc.

701 {
702 xc = xx0;
703 yc = yy0;
704 zc = z0;
705
706 if (m_w.empty()) return false;
707
708 double x0 = xx0;
709 double y0 = yy0;
710 double x1 = xx1;
711 double y1 = yy1;
712 if (m_polar) {
713 Cartesian2Internal(xx0, yy0, x0, y0);
714 Cartesian2Internal(xx1, yy1, x1, y1);
715 }
716 const double dx = x1 - x0;
717 const double dy = y1 - y0;
718 const double d2 = dx * dx + dy * dy;
719 // Check that the step length is non-zero.
720 if (d2 < Small) return false;
721 const double invd2 = 1. / d2;
722
723 // Check if a whole period has been crossed.
724 if ((m_perx && fabs(dx) >= m_sx) || (m_pery && fabs(dy) >= m_sy)) {
725 std::cerr << m_className << "::CrossedWire:\n"
726 << " Particle crossed more than one period.\n";
727 return false;
728 }
729
730 // Both coordinates are assumed to be located inside
731 // the drift area and inside a drift medium.
732 // This should have been checked before this call.
733
734 const double xm = 0.5 * (x0 + x1);
735 const double ym = 0.5 * (y0 + y1);
736 double dMin2 = 0.;
737 for (const auto& wire : m_w) {
738 double xw = wire.x;
739 if (m_perx) xw += m_sx * round((xm - xw) / m_sx);
740 double yw = wire.y;
741 if (m_pery) yw += m_sy * round((ym - yw) / m_sy);
742 // Calculate the smallest distance between track and wire.
743 const double xIn0 = dx * (xw - x0) + dy * (yw - y0);
744 // Check if the minimum is located before (x0, y0).
745 if (xIn0 < 0.) continue;
746 const double xIn1 = -(dx * (xw - x1) + dy * (yw - y1));
747 // Check if the minimum is located behind (x1, y1).
748 if (xIn1 < 0.) continue;
749 // Minimum is located between (x0, y0) and (x1, y1).
750 const double xw0 = xw - x0;
751 const double xw1 = xw - x1;
752 const double yw0 = yw - y0;
753 const double yw1 = yw - y1;
754 const double dw02 = xw0 * xw0 + yw0 * yw0;
755 const double dw12 = xw1 * xw1 + yw1 * yw1;
756 if (xIn1 * xIn1 * dw02 > xIn0 * xIn0 * dw12) {
757 dMin2 = dw02 - xIn0 * xIn0 * invd2;
758 } else {
759 dMin2 = dw12 - xIn1 * xIn1 * invd2;
760 }
761 // Add in the times nTrap to account for the trap radius.
762 const double r2 = wire.r * wire.r;
763 if (dMin2 < r2) {
764 // Wire has been crossed.
765 if (centre) {
766 if (m_polar) {
767 Internal2Cartesian(xw, yw, xc, yc);
768 } else {
769 xc = xw;
770 yc = yw;
771 }
772 } else {
773 // Find the point of intersection.
774 const double p = -xIn0 * invd2;
775 const double q = (dw02 - r2) * invd2;
776 const double s = sqrt(p * p - q);
777 const double t = std::min(-p + s, -p - s);
778 if (m_polar) {
779 Internal2Cartesian(x0 + t * dx, y0 + t * dy, xc, yc);
780 } else {
781 xc = x0 + t * dx;
782 yc = y0 + t * dy;
783 }
784 zc = z0 + t * (z1 - z0);
785 }
786 rc = wire.r;
787 if (m_polar) rc *= exp(wire.x);
788 return true;
789 }
790 }
791 return false;
792}
Heed::exp
DoubleAc exp(const DoubleAc &f)
Definition DoubleAc.cpp:377
Heed::sqrt
DoubleAc sqrt(const DoubleAc &f)
Definition DoubleAc.cpp:314

◆ ElectricField() [1/3]

std::array< double, 3 > Garfield::Component::ElectricField ( const double x,
const double y,
const double z )

Calculate the drift field [V/cm] at (x, y, z).

Definition at line 55 of file Component.cc.

43 {
44 double ex = 0., ey = 0., ez = 0.;
45 Medium* medium = nullptr;
46 int status = 0;
47 ElectricField(x, y, z, ex, ey, ez, medium, status);
48 std::array<double, 3> efield = {ex, ey, ez};
49 return efield;
50}
Garfield::ComponentAnalyticField::ElectricField
void ElectricField(const double x, const double y, const double z, double &ex, double &ey, double &ez, Medium *&m, int &status) override
Definition ComponentAnalyticField.hh:325

◆ ElectricField() [2/3]

void Garfield::ComponentAnalyticField::ElectricField ( const double x,
const double y,
const double z,
double & ex,
double & ey,
double & ez,
double & v,
Medium *& m,
int & status )
inlineoverridevirtual

Calculate the drift field [V/cm] and potential [V] at (x, y, z).

Implements Garfield::Component.

Definition at line 342 of file ComponentAnalyticField.hh.

344 {
345 m = nullptr;
346 // Calculate the field.
347 status = Field(x, y, z, ex, ey, ez, v, true);
348 // If the field is ok, get the medium.
349 if (status == 0) {
350 m = m_geometry ? m_geometry->GetMedium(x, y, z) : m_medium;
351 if (!m) {
352 status = -6;
353 } else if (!m->IsDriftable()) {
354 status = -5;
355 }
356 }
357 }
Garfield::Component::m_geometry
Geometry * m_geometry
Pointer to the geometry.
Definition Component.hh:362

◆ ElectricField() [3/3]

void Garfield::ComponentAnalyticField::ElectricField ( const double x,
const double y,
const double z,
double & ex,
double & ey,
double & ez,
Medium *& m,
int & status )
inlineoverridevirtual

Calculate the drift field at given point.

Parameters
x,y,zcoordinates [cm].
ex,ey,ezcomponents of the electric field [V/cm].
mpointer to the medium at this location.
statusstatus flag

Status flags: 

        0: Inside an active medium
      > 0: Inside a wire of type X
-4 ... -1: On the side of a plane where no wires are
       -5: Inside the mesh but not in an active medium
       -6: Outside the mesh
      -10: Unknown potential type (should not occur)
    other: Other cases (should not occur)

Implements Garfield::Component.

Definition at line 325 of file ComponentAnalyticField.hh.

326 {
327 m = nullptr;
328 // Calculate the field.
329 double v = 0.;
330 status = Field(x, y, z, ex, ey, ez, v, false);
331 // If the field is ok, get the medium.
332 if (status == 0) {
333 m = m_geometry ? m_geometry->GetMedium(x, y, z) : m_medium;
334 if (!m) {
335 status = -6;
336 } else if (!m->IsDriftable()) {
337 status = -5;
338 }
339 }
340 }

◆ ElectricFieldAtWire()

bool Garfield::ComponentAnalyticField::ElectricFieldAtWire ( const unsigned int iw,
double & ex,
double & ey )

Calculate the electric field at a given wire position, as if the wire itself were not there, but with the presence of its mirror images.

Definition at line 1933 of file ComponentAnalyticField.cc.

1934 {
1935 //-----------------------------------------------------------------------
1936 // FFIELD - Subroutine calculating the electric field at a given wire
1937 // position, as if the wire itself were not there but with
1938 // the presence of its mirror images.
1939 // VARIABLES : IW : wire number
1940 // EX, EY : x- and y-component of the electric field.
1941 // (Last changed on 27/ 1/96.)
1942 //-----------------------------------------------------------------------
1943 ex = ey = 0.;
1944 // Check the wire number.
1945 if (iw >= m_nWires) {
1946 std::cerr << m_className << "::ElectricFieldAtWire: Index out of range.\n";
1947 return false;
1948 }
1949 // Set the flags appropriately.
1950 std::vector<bool> cnalso(m_nWires, true);
1951 cnalso[iw] = false;
1952
1953 const double xpos = m_w[iw].x;
1954 const double ypos = m_w[iw].y;
1955 // Call the appropriate function.
1956 switch (m_cellType) {
1957 case A00:
1958 FieldAtWireA00(xpos, ypos, ex, ey, cnalso);
1959 break;
1960 case B1X:
1961 FieldAtWireB1X(xpos, ypos, ex, ey, cnalso);
1962 break;
1963 case B1Y:
1964 FieldAtWireB1Y(xpos, ypos, ex, ey, cnalso);
1965 break;
1966 case B2X:
1967 FieldAtWireB2X(xpos, ypos, ex, ey, cnalso);
1968 break;
1969 case B2Y:
1970 FieldAtWireB2Y(xpos, ypos, ex, ey, cnalso);
1971 break;
1972 case C10:
1973 FieldAtWireC10(xpos, ypos, ex, ey, cnalso);
1974 break;
1975 case C2X:
1976 FieldAtWireC2X(xpos, ypos, ex, ey, cnalso);
1977 break;
1978 case C2Y:
1979 FieldAtWireC2Y(xpos, ypos, ex, ey, cnalso);
1980 break;
1981 case C30:
1982 FieldAtWireC30(xpos, ypos, ex, ey, cnalso);
1983 break;
1984 case D10:
1985 FieldAtWireD10(xpos, ypos, ex, ey, cnalso);
1986 break;
1987 case D20:
1988 FieldAtWireD20(xpos, ypos, ex, ey, cnalso);
1989 break;
1990 case D30:
1991 FieldAtWireD30(xpos, ypos, ex, ey, cnalso);
1992 break;
1993 default:
1994 std::cerr << m_className << "::ElectricFieldAtWire:\n"
1995 << " Unknown cell type (id " << m_cellType << ")\n";
1996 return false;
1997 }
1998 // Correct for the equipotential planes.
1999 ex -= m_corvta;
2000 ey -= m_corvtb;
2001 if (m_polar) {
2002 const double r = exp(xpos);
2003 const double er = ex / r;
2004 const double ep = ey / r;
2005 const double ct = cos(ypos);
2006 const double st = sin(ypos);
2007 ex = +ct * er - st * ep;
2008 ey = +st * er + ct * ep;
2009 }
2010 return true;
2011}

Referenced by ForcesOnWire(), and WireDisplacement().

◆ EnableChargeCheck()

void Garfield::ComponentAnalyticField::EnableChargeCheck ( const bool on = true)
inline

Check the quality of the capacitance matrix inversion (default: off).

Definition at line 182 of file ComponentAnalyticField.hh.

182{ m_chargeCheck = on; }

◆ EnableDipoleTerms()

void Garfield::ComponentAnalyticField::EnableDipoleTerms ( const bool on = true)

Request dipole terms be included for each of the wires (default: off).

Definition at line 1597 of file ComponentAnalyticField.cc.

1597 {
1598
1599 m_cellset = false;
1600 m_sigset = false;
1601 m_dipole = on;
1602}

◆ EnableElectrostaticForce()

void Garfield::ComponentAnalyticField::EnableElectrostaticForce ( const bool on = true)
inline

Include the electrostatic force in the sag computation or not.

Definition at line 290 of file ComponentAnalyticField.hh.

290 {
291 m_useElectrostaticForce = on;
292 }

◆ EnableExtrapolation()

void Garfield::ComponentAnalyticField::EnableExtrapolation ( const bool on = true)
inline

Switch on/off extrapolation of electrostatic forces beyond the scanning area (default: off).

Definition at line 281 of file ComponentAnalyticField.hh.

281{ m_extrapolateForces = on; }

◆ EnableGravity()

void Garfield::ComponentAnalyticField::EnableGravity ( const bool on = true)
inline

Include gravity in the sag computation or not.

Definition at line 288 of file ComponentAnalyticField.hh.

288{ m_useGravitationalForce = on; }

◆ ForcesOnWire()

bool Garfield::ComponentAnalyticField::ForcesOnWire ( const unsigned int iw,
std::vector< double > & xMap,
std::vector< double > & yMap,
std::vector< std::vector< double > > & fxMap,
std::vector< std::vector< double > > & fyMap )

Calculate a table of the forces acting on a wire.

Parameters
iwindex of the wire
xMapcoordinates of the grid lines in x
yMapcoordinates of the grid lines in y
fxMapx-components of the force at the grid points
fyMapy-components of the force at the grid points

Definition at line 2077 of file ComponentAnalyticField.cc.

2080 {
2081 if (!m_cellset && !Prepare()) {
2082 std::cerr << m_className << "::ForcesOnWire: Cell not set up.\n";
2083 return false;
2084 }
2085
2086 if (iw >= m_nWires) {
2087 std::cerr << m_className << "::ForcesOnWire: Wire index out of range.\n";
2088 return false;
2089 }
2090 if (m_polar) {
2091 std::cerr << m_className << "::ForcesOnWire: Cannot handle polar cells.\n";
2092 return false;
2093 }
2094 const auto& wire = m_w[iw];
2095 // Compute a 'safe box' around the wire.
2096 double bxmin = m_perx ? wire.x - 0.5 * m_sx : 2 * m_xmin - m_xmax;
2097 double bxmax = m_perx ? wire.x + 0.5 * m_sx : 2 * m_xmax - m_xmin;
2098 double bymin = m_pery ? wire.y - 0.5 * m_sy : 2 * m_ymin - m_ymax;
2099 double bymax = m_pery ? wire.y + 0.5 * m_sy : 2 * m_ymax - m_ymin;
2100
2101 // If the initial area is almost zero in 1 direction, make it square.
2102 if (std::abs(bxmax - bxmin) < 0.1 * std::abs(bymax - bymin)) {
2103 bxmin = wire.x - 0.5 * std::abs(bymax - bymin);
2104 bxmax = wire.x + 0.5 * std::abs(bymax - bymin);
2105 } else if (std::abs(bymax - bymin) < 0.1 * std::abs(bxmax - bxmin)) {
2106 bymin = wire.y - 0.5 * std::abs(bxmax - bxmin);
2107 bymax = wire.y + 0.5 * std::abs(bxmax - bxmin);
2108 }
2109 const double dw = 2 * wire.r;
2110 // Scan the other wires.
2111 for (unsigned int j = 0; j < m_nWires; ++j) {
2112 if (j == iw) continue;
2113 const double xj = m_w[j].x;
2114 const double yj = m_w[j].y;
2115 const double dj = 2 * m_w[j].r;
2116 double xnear = m_perx ? xj - m_sx * round((xj - wire.x) / m_sx) : xj;
2117 double ynear = m_pery ? yj - m_sy * round((yj - wire.y) / m_sy) : yj;
2118 if (std::abs(xnear - wire.x) > std::abs(ynear - wire.y)) {
2119 if (xnear < wire.x) {
2120 bxmin = std::max(bxmin, xnear + dj + dw);
2121 if (m_perx) bxmax = std::min(bxmax, xnear + m_sx - dj - dw);
2122 } else {
2123 bxmax = std::min(bxmax, xnear - dj - dw);
2124 if (m_perx) bxmin = std::max(bxmin, xnear - m_sx + dj + dw);
2125 }
2126 } else {
2127 if (ynear < wire.y) {
2128 bymin = std::max({bymin, ynear - dj - dw, ynear + dj + dw});
2129 if (m_pery) bymax = std::min(bymax, ynear + m_sy - dj - dw);
2130 } else {
2131 bymax = std::min({bymax, ynear - dj - dw, ynear + dj + dw});
2132 if (m_pery) bymin = std::max(bymin, ynear - m_sy + dj + dw);
2133 }
2134 }
2135 }
2136 // Scan the planes.
2137 if (m_ynplan[0]) bxmin = std::max(bxmin, m_coplan[0] + dw);
2138 if (m_ynplan[1]) bxmax = std::min(bxmax, m_coplan[1] - dw);
2139 if (m_ynplan[2]) bymin = std::max(bymin, m_coplan[2] + dw);
2140 if (m_ynplan[3]) bymax = std::min(bymax, m_coplan[3] - dw);
2141
2142 // If there is a tube, check all corners.
2143 if (m_tube) {
2144 const double d2 = m_cotube2 - dw * dw;
2145 if (d2 < Small) {
2146 std::cerr << m_className << "::ForcesOnWire:\n Diameter of wire " << iw
2147 << " is too large compared to the tube.\n";
2148 return false;
2149 }
2150
2151 double corr = sqrt((bxmin * bxmin + bymin * bymin) / d2);
2152 if (corr > 1.) {
2153 bxmin /= corr;
2154 bymin /= corr;
2155 }
2156 corr = sqrt((bxmin * bxmin + bymax * bymax) / d2);
2157 if (corr > 1.) {
2158 bxmin /= corr;
2159 bymax /= corr;
2160 }
2161 corr = sqrt((bxmax * bxmax + bymin * bymin) / d2);
2162 if (corr > 1.) {
2163 bxmax /= corr;
2164 bymin /= corr;
2165 }
2166 corr = sqrt((bxmax * bxmax + bymax * bymax) / d2);
2167 if (corr > 1) {
2168 bxmax /= corr;
2169 bymax /= corr;
2170 }
2171 }
2172 // Make sure we found a reasonable 'safe area'.
2173 if ((bxmin - wire.x) * (wire.x - bxmax) <= 0 ||
2174 (bymin - wire.y) * (wire.y - bymax) <= 0) {
2175 std::cerr << m_className << "::ForcesOnWire:\n Unable to find an area "
2176 << "free of elements around wire " << iw << ".\n";
2177 return false;
2178 }
2179 // Now set a reasonable scanning range.
2180 double sxmin = bxmin;
2181 double sxmax = bxmax;
2182 double symin = bymin;
2183 double symax = bymax;
2184 if (m_scanRange == ScanningRange::User) {
2185 // User-specified range.
2186 sxmin = wire.x + m_xScanMin;
2187 symin = wire.y + m_yScanMin;
2188 sxmax = wire.x + m_xScanMax;
2189 symax = wire.y + m_yScanMax;
2190 } else if (m_scanRange == ScanningRange::FirstOrder) {
2191 // Get the field and force at the nominal position.
2192 double ex = 0., ey = 0.;
2193 ElectricFieldAtWire(iw, ex, ey);
2194 double fx = -ex * wire.e * Internal2Newton;
2195 double fy = -ey * wire.e * Internal2Newton;
2196 if (m_useGravitationalForce) {
2197 // Mass per unit length [kg / cm].
2198 const double m = 1.e-3 * wire.density * Pi * wire.r * wire.r;
2199 fx -= m_down[0] * GravitationalAcceleration * m;
2200 fy -= m_down[1] * GravitationalAcceleration * m;
2201 }
2202 const double u2 = wire.u * wire.u;
2203 const double shiftx =
2204 -125 * fx * u2 / (GravitationalAcceleration * wire.tension);
2205 const double shifty =
2206 -125 * fy * u2 / (GravitationalAcceleration * wire.tension);
2207 // If 0th order estimate of shift is not small.
2208 const double tol = 0.1 * wire.r;
2209 if (std::abs(shiftx) > tol || std::abs(shifty) > tol) {
2210 sxmin = std::max(bxmin, std::min(wire.x + m_scaleRange * shiftx,
2211 wire.x - shiftx / m_scaleRange));
2212 symin = std::max(bymin, std::min(wire.y + m_scaleRange * shifty,
2213 wire.y - shifty / m_scaleRange));
2214 sxmax = std::min(bxmax, std::max(wire.x + m_scaleRange * shiftx,
2215 wire.x - shiftx / m_scaleRange));
2216 symax = std::min(bymax, std::max(wire.y + m_scaleRange * shifty,
2217 wire.y - shifty / m_scaleRange));
2218 // If one is very small, make the area square within bounds.
2219 if (std::abs(sxmax - sxmin) < 0.1 * std::abs(symax - symin)) {
2220 sxmin = std::max(bxmin, wire.x - 0.5 * std::abs(symax - symin));
2221 sxmax = std::min(bxmax, wire.x + 0.5 * std::abs(symax - symin));
2222 } else if (std::abs(symax - symin) < 0.1 * std::abs(sxmax - sxmin)) {
2223 symin = std::max(bymin, wire.y - 0.5 * std::abs(sxmax - sxmin));
2224 symax = std::min(bymax, wire.y + 0.5 * std::abs(sxmax - sxmin));
2225 }
2226 }
2227 }
2228 if (m_debug) {
2229 std::cout << m_className << "::ForcesOnWire:\n";
2230 std::printf(" Free area %12.5e < x < %12.5e\n", bxmin, bxmax);
2231 std::printf(" %12.5e < y < %12.5e\n", bymin, bymax);
2232 std::printf(" Scan area %12.5e < x < %12.5e\n", sxmin, sxmax);
2233 std::printf(" %12.5e < y < %12.5e\n", symin, symax);
2234 }
2235
2236 xMap.resize(m_nScanX);
2237 const double stepx = (sxmax - sxmin) / (m_nScanX - 1);
2238 for (unsigned int i = 0; i < m_nScanX; ++i) {
2239 xMap[i] = sxmin + i * stepx;
2240 }
2241 yMap.resize(m_nScanY);
2242 const double stepy = (symax - symin) / (m_nScanY - 1);
2243 for (unsigned int i = 0; i < m_nScanY; ++i) {
2244 yMap[i] = symin + i * stepy;
2245 }
2246 // Save the original coordinates of the wire.
2247 const double x0 = wire.x;
2248 const double y0 = wire.y;
2249 // Prepare interpolation tables.
2250 fxMap.assign(m_nScanX, std::vector<double>(m_nScanY, 0.));
2251 fyMap.assign(m_nScanX, std::vector<double>(m_nScanY, 0.));
2252 bool ok = true;
2253 for (unsigned int i = 0; i < m_nScanX; ++i) {
2254 for (unsigned int j = 0; j < m_nScanY; ++j) {
2255 // Get the wire position for this shift.
2256 m_w[iw].x = xMap[i];
2257 m_w[iw].y = yMap[j];
2258 // Verify the current situation.
2259 if (!WireCheck()) {
2260 std::cerr << m_className << "::ForcesOnWire: Wire " << iw << ".\n"
2261 << " Scan involves a disallowed wire position.\n";
2262 ok = false;
2263 continue;
2264 }
2265 // Recompute the charges for this configuration.
2266 if (!Setup()) {
2267 std::cerr << m_className << "::ForcesOnWire: Wire " << iw << ".\n"
2268 << " Failed to compute charges at a scan point.\n";
2269 ok = false;
2270 continue;
2271 }
2272 // Compute the forces.
2273 double ex = 0., ey = 0.;
2274 ElectricFieldAtWire(iw, ex, ey);
2275 fxMap[i][j] = -ex * wire.e * Internal2Newton;
2276 fyMap[i][j] = -ey * wire.e * Internal2Newton;
2277 }
2278 }
2279 // Place the wire back in its shifted position.
2280 m_w[iw].x = x0;
2281 m_w[iw].y = y0;
2282 Setup();
2283 return ok;
2284}
Garfield::ComponentAnalyticField::ElectricFieldAtWire
bool ElectricFieldAtWire(const unsigned int iw, double &ex, double &ey)
Definition ComponentAnalyticField.cc:1933
Garfield::Component::m_debug
bool m_debug
Switch on/off debugging messages.
Definition Component.hh:371

Referenced by WireDisplacement().

◆ GetBoundingBox()

bool Garfield::ComponentAnalyticField::GetBoundingBox ( double & xmin,
double & ymin,
double & zmin,
double & xmax,
double & ymax,
double & zmax )
overridevirtual

Get the bounding box coordinates.

Reimplemented from Garfield::Component.

Definition at line 348 of file ComponentAnalyticField.cc.

350 {
351 // If a geometry is present, try to get the bounding box from there.
352 if (m_geometry) {
353 if (m_geometry->GetBoundingBox(x0, y0, z0, x1, y1, z1)) return true;
354 }
355 // Otherwise, return the cell dimensions.
356 return GetElementaryCell(x0, y0, z0, x1, y1, z1);
357}
Garfield::ComponentAnalyticField::GetElementaryCell
bool GetElementaryCell(double &x0, double &y0, double &z0, double &x1, double &y1, double &z1) override
Get the coordinates of the elementary cell.
Definition ComponentAnalyticField.cc:359

◆ GetCellType()

std::string Garfield::ComponentAnalyticField::GetCellType ( )
inline

Return the cell type. Cells are classified according to the number and orientation of planes, the presence of periodicities and the location of the wires as one of the following types:

A non-periodic cells with at most 1 x- and 1 y-plane B1X x-periodic cells without x-planes and at most 1 y-plane B1Y y-periodic cells without y-planes and at most 1 x-plane B2X cells with 2 x-planes and at most 1 y-plane B2Y cells with 2 y-planes and at most 1 x-plane C1 doubly periodic cells without planes C2X doubly periodic cells with x-planes C2Y doubly periodic cells with y-planes C3 double periodic cells with x- and y-planes D1 round tubes without axial periodicity D2 round tubes with axial periodicity D3 polygonal tubes without axial periodicity

Definition at line 152 of file ComponentAnalyticField.hh.

152 {
153 if (!m_cellset) {
154 if (CellCheck()) CellType();
155 }
156 return GetCellType(m_cellType);
157 }

Referenced by GetCellType(), and PrintCell().

◆ GetElementaryCell()

bool Garfield::ComponentAnalyticField::GetElementaryCell ( double & xmin,
double & ymin,
double & zmin,
double & xmax,
double & ymax,
double & zmax )
overridevirtual

Get the coordinates of the elementary cell.

Reimplemented from Garfield::Component.

Definition at line 359 of file ComponentAnalyticField.cc.

361 {
362 if (!m_cellset && !Prepare()) return false;
363 if (m_polar) {
364 double rmax, thetamax;
365 Internal2Polar(m_xmax, m_ymax, rmax, thetamax);
366 x0 = -rmax;
367 y0 = -rmax;
368 x1 = +rmax;
369 y1 = +rmax;
370 } else {
371 x0 = m_xmin;
372 y0 = m_ymin;
373 x1 = m_xmax;
374 y1 = m_ymax;
375 }
376 z0 = m_zmin;
377 z1 = m_zmax;
378 return true;
379}

Referenced by GetBoundingBox().

◆ GetGravity()

void Garfield::ComponentAnalyticField::GetGravity ( double & dx,
double & dy,
double & dz ) const

Get the gravity orientation.

Definition at line 2069 of file ComponentAnalyticField.cc.

2070 {
2071
2072 dx = m_down[0];
2073 dy = m_down[1];
2074 dz = m_down[2];
2075}

◆ GetMedium()

Medium * Garfield::ComponentAnalyticField::GetMedium ( const double x,
const double y,
const double z )
overridevirtual

Get the medium at a given location (x, y, z).

Reimplemented from Garfield::Component.

Definition at line 283 of file ComponentAnalyticField.cc.

284 {
285 if (m_geometry) return m_geometry->GetMedium(xin, yin, zin);
286
287 // Make sure the cell is prepared.
288 if (!m_cellset && !Prepare()) return nullptr;
289
290 double xpos = xin, ypos = yin;
291 if (m_polar) Cartesian2Internal(xin, yin, xpos, ypos);
292 // In case of periodicity, move the point into the basic cell.
293 if (m_perx) xpos -= m_sx * round(xin / m_sx);
294
295 double arot = 0.;
296 if (m_pery && m_tube) {
297 Cartesian2Polar(xin, yin, xpos, ypos);
298 arot = RadToDegree * m_sy * round(DegreeToRad * ypos / m_sy);
299 ypos -= arot;
300 Polar2Cartesian(xpos, ypos, xpos, ypos);
301 } else if (m_pery) {
302 ypos -= m_sy * round(ypos / m_sy);
303 }
304
305 // Move the point to the correct side of the plane.
306 if (m_perx && m_ynplan[0] && xpos <= m_coplan[0]) xpos += m_sx;
307 if (m_perx && m_ynplan[1] && xpos >= m_coplan[1]) xpos -= m_sx;
308 if (m_pery && m_ynplan[2] && ypos <= m_coplan[2]) ypos += m_sy;
309 if (m_pery && m_ynplan[3] && ypos >= m_coplan[3]) ypos -= m_sy;
310
311 // In case (xpos, ypos) is located behind a plane there is no field.
312 if (m_tube) {
313 if (!InTube(xpos, ypos, m_cotube, m_ntube)) return nullptr;
314 } else {
315 if ((m_ynplan[0] && xpos < m_coplan[0]) ||
316 (m_ynplan[1] && xpos > m_coplan[1]) ||
317 (m_ynplan[2] && ypos < m_coplan[2]) ||
318 (m_ynplan[3] && ypos > m_coplan[3])) {
319 return nullptr;
320 }
321 }
322
323 // If (xpos, ypos) is within a wire, there is no field either.
324 for (const auto& wire : m_w) {
325 double dx = xpos - wire.x;
326 double dy = ypos - wire.y;
327 // Correct for periodicities.
328 if (m_perx) dx -= m_sx * round(dx / m_sx);
329 if (m_pery) dy -= m_sy * round(dy / m_sy);
330 // Check the actual position.
331 if (dx * dx + dy * dy < wire.r * wire.r) return nullptr;
332 }
333 return m_medium;
334}

◆ GetNumberOfPlanesPhi()

unsigned int Garfield::ComponentAnalyticField::GetNumberOfPlanesPhi ( ) const

Get the number of equipotential planes at constant phi.

Definition at line 1828 of file ComponentAnalyticField.cc.

1828 {
1829 if (!m_polar) {
1830 return 0;
1831 } else if (m_ynplan[2] && m_ynplan[3]) {
1832 return 2;
1833 } else if (m_ynplan[2] || m_ynplan[3]) {
1834 return 1;
1835 }
1836 return 0;
1837}

◆ GetNumberOfPlanesR()

unsigned int Garfield::ComponentAnalyticField::GetNumberOfPlanesR ( ) const

Get the number of equipotential planes at constant radius.

Definition at line 1817 of file ComponentAnalyticField.cc.

1817 {
1818 if (!m_polar) {
1819 return 0;
1820 } else if (m_ynplan[0] && m_ynplan[1]) {
1821 return 2;
1822 } else if (m_ynplan[0] || m_ynplan[1]) {
1823 return 1;
1824 }
1825 return 0;
1826}

◆ GetNumberOfPlanesX()

unsigned int Garfield::ComponentAnalyticField::GetNumberOfPlanesX ( ) const

Get the number of equipotential planes at constant x.

Definition at line 1795 of file ComponentAnalyticField.cc.

1795 {
1796 if (m_polar) {
1797 return 0;
1798 } else if (m_ynplan[0] && m_ynplan[1]) {
1799 return 2;
1800 } else if (m_ynplan[0] || m_ynplan[1]) {
1801 return 1;
1802 }
1803 return 0;
1804}

◆ GetNumberOfPlanesY()

unsigned int Garfield::ComponentAnalyticField::GetNumberOfPlanesY ( ) const

Get the number of equipotential planes at constant y.

Definition at line 1806 of file ComponentAnalyticField.cc.

1806 {
1807 if (m_polar) {
1808 return 0;
1809 } else if (m_ynplan[2] && m_ynplan[3]) {
1810 return 2;
1811 } else if (m_ynplan[2] || m_ynplan[3]) {
1812 return 1;
1813 }
1814 return 0;
1815}

◆ GetNumberOfWires()

unsigned int Garfield::ComponentAnalyticField::GetNumberOfWires ( ) const
inline

Get the number of wires.

Definition at line 185 of file ComponentAnalyticField.hh.

185{ return m_nWires; }

◆ GetPeriodicityPhi()

bool Garfield::ComponentAnalyticField::GetPeriodicityPhi ( double & s)

Get the periodicity [degree] in phi.

Definition at line 1676 of file ComponentAnalyticField.cc.

1676 {
1677 if (!m_periodic[1] || (!m_polar && !m_tube)) {
1678 s = 0.;
1679 return false;
1680 }
1681
1682 s = RadToDegree * m_sy;
1683 return true;
1684}
Garfield::Component::m_periodic
std::array< bool, 3 > m_periodic
Simple periodicity in x, y, z.
Definition Component.hh:374

◆ GetPeriodicityX()

bool Garfield::ComponentAnalyticField::GetPeriodicityX ( double & s)

Get the periodic length in the x-direction.

Definition at line 1656 of file ComponentAnalyticField.cc.

1656 {
1657 if (!m_periodic[0] || m_polar) {
1658 s = 0.;
1659 return false;
1660 }
1661
1662 s = m_sx;
1663 return true;
1664}

◆ GetPeriodicityY()

bool Garfield::ComponentAnalyticField::GetPeriodicityY ( double & s)

Get the periodic length in the y-direction.

Definition at line 1666 of file ComponentAnalyticField.cc.

1666 {
1667 if (!m_periodic[1] || m_polar) {
1668 s = 0.;
1669 return false;
1670 }
1671
1672 s = m_sy;
1673 return true;
1674}

◆ GetPlanePhi()

bool Garfield::ComponentAnalyticField::GetPlanePhi ( const unsigned int i,
double & phi,
double & voltage,
std::string & label ) const

Retrieve the parameters of a plane at constant phi.

Definition at line 1909 of file ComponentAnalyticField.cc.

1911 {
1912 if (!m_polar || i >= 2 || (i == 1 && !m_ynplan[3])) {
1913 std::cerr << m_className << "::GetPlanePhi: Index out of range.\n";
1914 return false;
1915 }
1916
1917 phi = RadToDegree * m_coplan[i + 2];
1918 voltage = m_vtplan[i + 2];
1919 label = m_planes[i + 2].type;
1920 return true;
1921}

◆ GetPlaneR()

bool Garfield::ComponentAnalyticField::GetPlaneR ( const unsigned int i,
double & r,
double & voltage,
std::string & label ) const

Retrieve the parameters of a plane at constant radius.

Definition at line 1895 of file ComponentAnalyticField.cc.

1897 {
1898 if (!m_polar || i >= 2 || (i == 1 && !m_ynplan[1])) {
1899 std::cerr << m_className << "::GetPlaneR: Index out of range.\n";
1900 return false;
1901 }
1902
1903 r = exp(m_coplan[i]);
1904 voltage = m_vtplan[i];
1905 label = m_planes[i].type;
1906 return true;
1907}

◆ GetPlaneX()

bool Garfield::ComponentAnalyticField::GetPlaneX ( const unsigned int i,
double & x,
double & voltage,
std::string & label ) const

Retrieve the parameters of a plane at constant x.

Definition at line 1867 of file ComponentAnalyticField.cc.

1869 {
1870 if (m_polar || i >= 2 || (i == 1 && !m_ynplan[1])) {
1871 std::cerr << m_className << "::GetPlaneX: Index out of range.\n";
1872 return false;
1873 }
1874
1875 x = m_coplan[i];
1876 voltage = m_vtplan[i];
1877 label = m_planes[i].type;
1878 return true;
1879}

◆ GetPlaneY()

bool Garfield::ComponentAnalyticField::GetPlaneY ( const unsigned int i,
double & y,
double & voltage,
std::string & label ) const

Retrieve the parameters of a plane at constant y.

Definition at line 1881 of file ComponentAnalyticField.cc.

1883 {
1884 if (m_polar || i >= 2 || (i == 1 && !m_ynplan[3])) {
1885 std::cerr << m_className << "::GetPlaneY: Index out of range.\n";
1886 return false;
1887 }
1888
1889 y = m_coplan[i + 2];
1890 voltage = m_vtplan[i + 2];
1891 label = m_planes[i + 2].type;
1892 return true;
1893}

◆ GetTube()

bool Garfield::ComponentAnalyticField::GetTube ( double & r,
double & voltage,
int & nEdges,
std::string & label ) const

Retrieve the tube parameters.

Definition at line 1923 of file ComponentAnalyticField.cc.

1924 {
1925 if (!m_tube) return false;
1926 r = m_cotube;
1927 voltage = m_vttube;
1928 nEdges = m_ntube;
1929 label = m_planes[4].type;
1930 return true;
1931}

◆ GetVoltageRange()

bool Garfield::ComponentAnalyticField::GetVoltageRange ( double & vmin,
double & vmax )
overridevirtual

Calculate the voltage range [V].

Implements Garfield::Component.

Definition at line 336 of file ComponentAnalyticField.cc.

336 {
337 // Make sure the cell is prepared.
338 if (!m_cellset && !Prepare()) {
339 std::cerr << m_className << "::GetVoltageRange: Cell not set up.\n";
340 return false;
341 }
342
343 pmin = m_vmin;
344 pmax = m_vmax;
345 return true;
346}

◆ GetWire()

bool Garfield::ComponentAnalyticField::GetWire ( const unsigned int i,
double & x,
double & y,
double & diameter,
double & voltage,
std::string & label,
double & length,
double & charge,
int & ntrap ) const

Retrieve the parameters of a wire.

Definition at line 1839 of file ComponentAnalyticField.cc.

1842 {
1843 if (i >= m_nWires) {
1844 std::cerr << m_className << "::GetWire: Index out of range.\n";
1845 return false;
1846 }
1847
1848 if (m_polar) {
1849 double r = 0., theta = 0.;
1850 Internal2Polar(m_w[i].x, m_w[i].y, r, theta);
1851 x = r;
1852 y = theta;
1853 diameter = 2 * m_w[i].r * r;
1854 } else {
1855 x = m_w[i].x;
1856 y = m_w[i].y;
1857 diameter = 2 * m_w[i].r;
1858 }
1859 voltage = m_w[i].v;
1860 label = m_w[i].type;
1861 length = m_w[i].u;
1862 charge = m_w[i].e;
1863 ntrap = m_w[i].nTrap;
1864 return true;
1865}

◆ InTrapRadius()

bool Garfield::ComponentAnalyticField::InTrapRadius ( const double q0,
const double x0,
const double y0,
const double z0,
double & xw,
double & yw,
double & rw )
overridevirtual

Determine whether a particle is inside the trap radius of a wire.

Parameters
q0charge of the particle [in elementary charges].
x0,y0,z0position [cm] of the particle.
xw,ywcoordinates of the wire (if applicable).
rwradius of the wire (if applicable).

Reimplemented from Garfield::Component.

Definition at line 794 of file ComponentAnalyticField.cc.

797 {
798
799 double x0 = xin;
800 double y0 = yin;
801 if (m_polar) {
802 Cartesian2Internal(xin, yin, x0, y0);
803 }
804 // In case of periodicity, move the point into the basic cell.
805 int nX = 0, nY = 0, nPhi = 0;
806 if (m_perx) {
807 nX = int(round(x0 / m_sx));
808 x0 -= m_sx * nX;
809 }
810 if (m_pery && m_tube) {
811 Cartesian2Polar(xin, yin, x0, y0);
812 nPhi = int(round(DegreeToRad * y0 / m_sy));
813 y0 -= RadToDegree * m_sy * nPhi;
814 Polar2Cartesian(x0, y0, x0, y0);
815 } else if (m_pery) {
816 nY = int(round(y0 / m_sy));
817 y0 -= m_sy * nY;
818 }
819 // Move the point to the correct side of the plane.
820 std::array<bool, 4> shift = {false, false, false, false};
821 if (m_perx && m_ynplan[0] && x0 <= m_coplan[0]) {
822 x0 += m_sx;
823 shift[0] = true;
824 }
825 if (m_perx && m_ynplan[1] && x0 >= m_coplan[1]) {
826 x0 -= m_sx;
827 shift[1] = true;
828 }
829 if (m_pery && m_ynplan[2] && y0 <= m_coplan[2]) {
830 y0 += m_sy;
831 shift[2] = true;
832 }
833 if (m_pery && m_ynplan[3] && y0 >= m_coplan[3]) {
834 y0 -= m_sy;
835 shift[3] = true;
836 }
837
838 double r2min = std::numeric_limits<double>::max();
839 bool inside = false;
840 for (const auto& wire : m_w) {
841 // Skip wires with the wrong charge.
842 if (qin * wire.e > 0.) continue;
843 const double dxw0 = wire.x - x0;
844 const double dyw0 = wire.y - y0;
845 const double r2 = dxw0 * dxw0 + dyw0 * dyw0;
846 const double rTrap = wire.r * wire.nTrap;
847 if (r2 > rTrap * rTrap || r2 > r2min) continue;
848 inside = true;
849 r2min = r2;
850 xw = wire.x;
851 yw = wire.y;
852 rw = wire.r;
853 if (shift[0]) xw -= m_sx;
854 if (shift[1]) xw += m_sx;
855 if (shift[2]) yw -= m_sy;
856 if (shift[3]) yw += m_sy;
857 if (m_pery && m_tube) {
858 double rhow, phiw;
859 Cartesian2Polar(xw, yw, rhow, phiw);
860 phiw += RadToDegree * m_sy * nPhi;
861 Polar2Cartesian(rhow, phiw, xw, yw);
862 } else if (m_pery) {
863 y0 += m_sy * nY;
864 }
865 if (m_perx) xw += m_sx * nX;
866 if (m_polar) {
867 Internal2Cartesian(xw, yw, xw, yw);
868 rw *= exp(wire.x);
869 }
870 if (m_debug) {
871 std::cout << m_className << "::InTrapRadius: (" << xin << ", "
872 << yin << ", " << zin << ")" << " within trap radius.\n";
873 }
874 }
875 return inside;
876}

◆ IsPolar()

bool Garfield::ComponentAnalyticField::IsPolar ( ) const
inline

Are polar coordinates being used?

Definition at line 118 of file ComponentAnalyticField.hh.

118{ return m_polar; }

◆ MultipoleMoments()

bool Garfield::ComponentAnalyticField::MultipoleMoments ( const unsigned int iw,
const unsigned int order = 4,
const bool print = false,
const bool plot = false,
const double rmult = 1.,
const double eps = 1.e-4,
const unsigned int nMaxIter = 20 )

Calculate multipole moments for a given wire.

Parameters
iwIndex of the wire.
orderOrder of the highest multipole moment.
printPrint information about the fitting process.
plotPlot the potential and multipole fit around the wire.
rmultDistance in multiples of the wire radius at which the decomposition is to be carried out.
epsUsed in the fit for calculating the covariance matrix.
nMaxIterMaximum number of iterations in the fit.

Definition at line 10606 of file ComponentAnalyticField.cc.

10609 {
10610
10611 //-----------------------------------------------------------------------
10612 // EFMWIR - Computes the dipole moment of a given wire.
10613 //-----------------------------------------------------------------------
10614 if (!m_cellset && !Prepare()) return false;
10615 // Check input parameters.
10616 if (iw >= m_nWires) {
10617 std::cerr << m_className << "::MultipoleMoments:\n"
10618 << " Wire index out of range.\n";
10619 return false;
10620 }
10621 if (eps <= 0.) {
10622 std::cerr << m_className << "::MultipoleMoments:\n"
10623 << " Epsilon must be positive.\n";
10624 return false;
10625 }
10626 if (nPoles == 0) {
10627 std::cerr << m_className << "::MultipoleMoments:\n"
10628 << " Multipole order out of range.\n";
10629 return false;
10630 }
10631 if (rmult <= 0.) {
10632 std::cerr << m_className << "::MultipoleMoments:\n"
10633 << " Radius multiplication factor out of range.\n";
10634 return false;
10635 }
10636
10637 const double xw = m_w[iw].x;
10638 const double yw = m_w[iw].y;
10639 // Set the radius of the wire to 0.
10640 const double rw = m_w[iw].r;
10641 m_w[iw].r = 0.;
10642
10643 // Loop around the wire.
10644 constexpr unsigned int nPoints = 20000;
10645 std::vector<double> angle(nPoints, 0.);
10646 std::vector<double> volt(nPoints, 0.);
10647 std::vector<double> weight(nPoints, 1.);
10648 for (unsigned int i = 0; i < nPoints; ++i) {
10649 // Set angle around wire.
10650 angle[i] = TwoPi * (i + 1.) / nPoints;
10651 // Compute E field, make sure the point is in a free region.
10652 const double x = xw + rmult * rw * cos(angle[i]);
10653 const double y = yw + rmult * rw * sin(angle[i]);
10654 double ex = 0., ey = 0., ez = 0., v = 0.;
10655 if (Field(x, y, 0., ex, ey, ez, v, true) != 0) {
10656 std::cerr << m_className << "::MultipoleMoments:\n"
10657 << " Unexpected location code. Computation stopped.\n";
10658 m_w[iw].r = rw;
10659 return false;
10660 }
10661 // Assign the result to the fitting array.
10662 volt[i] = v;
10663 }
10664 // Restore the wire diameter.
10665 m_w[iw].r = rw;
10666
10667 // Determine the maximum, minimum and average.
10668 double vmin = *std::min_element(volt.cbegin(), volt.cend());
10669 double vmax = *std::max_element(volt.cbegin(), volt.cend());
10670 double vave = std::accumulate(volt.cbegin(), volt.cend(), 0.) / nPoints;
10671 // Subtract the wire potential to centre the data more or less.
10672 for (unsigned int i = 0; i < nPoints; ++i) volt[i] -= vave;
10673 vmax -= vave;
10674 vmin -= vave;
10675
10676 // Perform the fit.
10677 const double vm = 0.5 * fabs(vmin) + fabs(vmax);
10678 double chi2 = 1.e-6 * nPoints * vm * vm;
10679 const double dist = 1.e-3 * (1. + vm);
10680 const unsigned int nPar = 2 * nPoles + 1;
10681 std::vector<double> pars(nPar, 0.);
10682 std::vector<double> epar(nPar, 0.);
10683 pars[0] = 0.5 * (vmax + vmin);
10684 for (unsigned int i = 1; i <= nPoles; ++i) {
10685 pars[2 * i - 1] = 0.5 * (vmax - vmin);
10686 pars[2 * i] = 0.;
10687 }
10688
10689 auto f = [nPoles](const double x, const std::vector<double>& par) {
10690 // EFMFUN
10691 // Sum the series, initial value is the monopole term.
10692 double sum = par[0];
10693 for (unsigned int k = 1; k <= nPoles; ++k) {
10694 // Obtain the Legendre polynomial of this order and add to the series.
10695 const float cphi = cos(x - par[2 * k]);
10696 sum += par[2 * k - 1] * sqrt(k + 0.5) * Numerics::Legendre(k, cphi);
10697 }
10698 return sum;
10699 };
10700
10701 if (!Numerics::LeastSquaresFit(f, pars, epar, angle, volt, weight,
10702 nMaxIter, dist, chi2, eps, m_debug, print)) {
10703 std::cerr << m_className << "::MultipoleMoments:\n"
10704 << " Fitting the multipoles failed; computation stopped.\n";
10705 }
10706 // Plot the result of the fit.
10707 if (plot) {
10708 const std::string name = ViewBase::FindUnusedCanvasName("cMultipole");
10709 TCanvas* cfit = new TCanvas(name.c_str(), "Multipoles");
10710 cfit->SetGridx();
10711 cfit->SetGridy();
10712 cfit->DrawFrame(0., vmin, TwoPi, vmax,
10713 ";Angle around the wire [rad]; Potential - average [V]");
10714 TGraph graph;
10715 graph.SetLineWidth(2);
10716 graph.SetLineColor(kBlack);
10717 graph.DrawGraph(angle.size(), angle.data(), volt.data(), "lsame");
10718 // Sum of contributions.
10719 constexpr unsigned int nP = 1000;
10720 std::array<double, nP> xp;
10721 std::array<double, nP> yp;
10722 for (unsigned int i = 0; i < nP; ++i) {
10723 xp[i] = TwoPi * (i + 1.) / nP;
10724 yp[i] = f(xp[i], pars);
10725 }
10726 graph.SetLineColor(kViolet + 3);
10727 graph.DrawGraph(nP, xp.data(), yp.data(), "lsame");
10728 // Individual contributions.
10729 std::vector<double> parres = pars;
10730 for (unsigned int i = 1; i <= nPoles; ++i) parres[2 * i - 1] = 0.;
10731 for (unsigned int j = 1; j <= nPoles; ++j) {
10732 parres[2 * j - 1] = pars[2 * j - 1];
10733 for (unsigned int i = 0; i < nP; ++i) {
10734 yp[i] = f(xp[i], parres);
10735 }
10736 parres[2 * j - 1] = 0.;
10737 graph.SetLineColor(kAzure + j);
10738 graph.DrawGraph(nP, xp.data(), yp.data(), "lsame");
10739 }
10740 gPad->Update();
10741 }
10742
10743 // Print the results.
10744 std::cout << m_className << "::MultipoleMoments:\n"
10745 << " Multipole moments for wire " << iw << ":\n"
10746 << " Moment Value Angle [degree]\n";
10747 std::printf(" %6u %15.8f Arbitrary\n", 0, vave);
10748 for (unsigned int i = 1; i <= nPoles; ++i) {
10749 // Remove radial term from the multipole moment.
10750 const double val = pow(rmult * rw, i) * pars[2 * i - 1];
10751 const double phi = RadToDegree * fmod(pars[2 * i], Pi);
10752 std::printf(" %6u %15.8f %15.8f\n", i, val, phi);
10753 }
10754 return true;
10755}
Garfield::ViewBase::FindUnusedCanvasName
static std::string FindUnusedCanvasName(const std::string &s)
Find an unused canvas name.
Definition ViewBase.cc:337
Garfield::Numerics::Legendre
double Legendre(const unsigned int n, const double x)
Legendre polynomials.
Definition Numerics.hh:162
Garfield::Numerics::LeastSquaresFit
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.
Definition Numerics.cc:1888
Heed::pow
DoubleAc pow(const DoubleAc &f, double p)
Definition DoubleAc.cpp:337

◆ OptimiseOnGrid()

bool Garfield::ComponentAnalyticField::OptimiseOnGrid ( const std::vector< std::string > & groups,
const std::string & field_function,
const double target,
const double x0,
const double y0,
const double x1,
const double y1,
const unsigned int nX = 10,
const unsigned int nY = 10,
const bool print = true )

Vary the potential of selected electrodes to match (a function of) the potential or field on an x-y (or r-phi) grid of points.

Definition at line 10936 of file ComponentAnalyticField.cc.

10941 {
10942 // -----------------------------------------------------------------------
10943 // OPTSET - Routine attempting to find proper voltage settings.
10944 // -----------------------------------------------------------------------
10945
10946 if (nX * nY < 1) {
10947 std::cerr << m_className << "::OptimiseOnGrid: Number of points < 1.\n";
10948 return false;
10949 }
10950
10951 if (groups.empty()) {
10952 std::cerr << m_className << "::OptimiseOnGrid: No electrode groups.\n";
10953 return false;
10954 }
10955
10956 // Use a constant target value for now.
10957 std::string target_function = std::to_string(target);
10958
10959 if (m_debug) {
10960 std::cout << m_className << "::OptimiseOnGrid:\n"
10961 << " The function " << field_function << "\n"
10962 << " has to approximate the value of the function "
10963 << target_function << ".\n";
10964 }
10965
10966 // Convert the field function.
10967 std::vector<std::string> variables;
10968 if (m_polar) {
10969 variables = {"r", "phi", "er", "ephi", "e", "v"};
10970 } else {
10971 variables = {"x", "y", "ex", "ey", "e", "v"};
10972 }
10973 auto fPtr = MakeFunction("fField", field_function, variables);
10974 if (!fPtr) {
10975 std::cerr << m_className << "::OptimiseOnGrid: "
10976 << "Error in the field function.\n";
10977 return false;
10978 }
10979 typedef std::function<double(const std::vector<double>&)> Fcn;
10980 Fcn& fField = *((Fcn *) fPtr);
10981
10982 // TODO: add option to use the current average of the field function
10983 // as target value.
10984
10985 // Convert the target function.
10986 if (m_polar) {
10987 variables = {"r", "phi"};
10988 } else {
10989 variables = {"x", "y"};
10990 }
10991 fPtr = MakeFunction("fTarget", target_function, variables);
10992 if (!fPtr) {
10993 std::cerr << m_className << "::OptimiseOnGrid: "
10994 << "Error in the target function.\n";
10995 return false;
10996 }
10997 Fcn& fTarget = *((Fcn *) fPtr);
10998
10999 const size_t nP = nX * nY;
11000 std::vector<double> xFit(nP);
11001 std::iota(xFit.begin(), xFit.end(), 0);
11002 std::vector<double> yFit(nP);
11003 std::vector<double> wFit(nP, 1.);
11004 const double dx = m_polar ? (exp(x1) - exp(x0)) / double(nX - 1) :
11005 (x1 - x0) / double(nX - 1);
11006 const double dy = (y1 - y0) / double(nY - 1);
11007 for (unsigned int i = 0; i < nX; ++i) {
11008 for (unsigned int j = 0; j < nY; ++j) {
11009 const unsigned int k = i + nX * j;
11010 // Grid position.
11011 std::vector<double> var(2, 0.);
11012 var[0] = m_polar ? log(exp(x0) + i * dx) : x0 + i * dx;
11013 var[1] = y0 + j * dy;
11014 if (m_polar) Internal2Polar(var[0], var[1], var[0], var[1]);
11015 // Evaluate the position dependent target function.
11016 yFit[k] = fTarget(var);
11017 // Evaluate the position dependent weighting function.
11018 // TODO
11019 }
11020 }
11021
11022 // Initialise the fit parameters.
11023 std::vector<double> vw0;
11024 std::array<double, 5> vp0;
11025 std::vector<double> aFit;
11026 std::vector<std::vector<unsigned int> > wiresInGroup;
11027 std::vector<std::vector<unsigned int> > planesInGroup;
11028 InitialiseFitParameters(groups, vw0, vp0, aFit, wiresInGroup, planesInGroup);
11029 if (aFit.empty()) {
11030 std::cerr << m_className << "::OptimiseOnGrid: "
11031 << "Setting fitting parameters failed.\n";
11032 return false;
11033 }
11034
11035 auto fFit = [&](const double s, const std::vector<double>& par) {
11036 // First set the potentials.
11037 std::vector<double> vw = vw0;
11038 std::array<double, 5> vp = vp0;
11039 const size_t nPar = par.size();
11040 for (size_t i = 0; i < nPar; ++i) {
11041 for (unsigned int iw : wiresInGroup[i]) vw[iw] += par[i];
11042 for (size_t ip : planesInGroup[i]) vp[ip] += par[i];
11043 }
11044 // Next reconstruct the charges.
11045 if (!Update(vw, vp)) return 0.;
11046
11047 std::vector<double> var(6, 0.);
11048 const int j = int(s) / nX;
11049 const int i = int(s) - nX * j;
11050 var[0] = m_polar ? log(exp(x0) + i * dx) : x0 + i * dx;
11051 var[1] = y0 + j * dy;
11052 double ez = 0.;
11053 Field(var[0], var[1], 0., var[2], var[3], ez, var[5], true);
11054 var[4] = sqrt(var[2] * var[2] + var[3] * var[3]);
11055 // TODO: drift time, diffusion, gain, ...
11056 // Transform to polar if needed.
11057 if (m_polar) {
11058 Internal2Polar(var[0], var[1], var[0], var[1]);
11059 var[2] = var[2] / var[0];
11060 var[3] = var[3] / var[0];
11061 var[4] = var[4] / var[0];
11062 }
11063 // Calculate the field function with this field.
11064 return fField(var);
11065 };
11066
11067 double chi2 = 0.;
11068 std::vector<double> eFit(aFit.size(), 0.);
11069 // Carry out the fitting itself.
11070 if (!Numerics::LeastSquaresFit(fFit, aFit, eFit, xFit, yFit, wFit,
11071 m_optNitmax, m_optDist, chi2, m_optEps,
11072 print, m_debug)) {
11073 std::cerr << m_className << "::OptimiseOnGrid:\n"
11074 << " The new potentials do not fulfill your requirements.\n";
11075 return false;
11076 }
11077
11078 // Calculate the charges for the final result.
11079 std::vector<double> vw = vw0;
11080 std::array<double, 5> vp = vp0;
11081 for (size_t i = 0; i < aFit.size(); ++i) {
11082 for (unsigned int iw : wiresInGroup[i]) vw[iw] += aFit[i];
11083 for (size_t ip : planesInGroup[i]) vp[ip] += aFit[i];
11084 }
11085 if (!Update(vw, vp)) {
11086 std::cerr << m_className << "::OptimiseOnGrid:\n"
11087 << " Failed to compute the charges for the final settings.\n";
11088 return false;
11089 }
11090 return true;
11091}

◆ OptimiseOnTrack()

bool Garfield::ComponentAnalyticField::OptimiseOnTrack ( const std::vector< std::string > & groups,
const std::string & field_function,
const double target,
const double x0,
const double y0,
const double x1,
const double y1,
const unsigned int nP = 20,
const bool print = true )

Vary the potential of selected electrodes to match (a function of) the potential or field along a given line to a target.

Parameters
groupsIdentifier of the electrodes for which to vary the potential.
field_functionA function of the coordinates (x, y or r, phi), the electrostatic field (ex, ey, e) and potential (v).
targetValue of the field function to be reproduced.
x0,y0Starting point of the line.
x1,y1End point of the line.
nPNumber of points on the line.
printFlag to print out information during each fit cycle.

Definition at line 10780 of file ComponentAnalyticField.cc.

10785 {
10786
10787 // -----------------------------------------------------------------------
10788 // OPTSET - Routine attempting to find proper voltage settings.
10789 // -----------------------------------------------------------------------
10790 if (nP < 1) {
10791 std::cerr << m_className << "::OptimiseOnTrack: Number of points < 1.\n";
10792 return false;
10793 }
10794
10795 if (groups.empty()) {
10796 std::cerr << m_className << "::OptimiseOnTrack: No electrode groups.\n";
10797 return false;
10798 }
10799
10800 // Use a constant target value for now.
10801 std::string target_function = std::to_string(target);
10802
10803 if (m_debug) {
10804 std::cout << m_className << "::OptimiseOnTrack:\n"
10805 << " The function " << field_function << "\n"
10806 << " has to approximate the value of the function "
10807 << target_function << ".\n";
10808 }
10809
10810 // Convert the field function.
10811 std::vector<std::string> variables;
10812 if (m_polar) {
10813 variables = {"r", "phi", "er", "ephi", "e", "v"};
10814 } else {
10815 variables = {"x", "y", "ex", "ey", "e", "v"};
10816 }
10817 auto fPtr = MakeFunction("fField", field_function, variables);
10818 if (!fPtr) {
10819 std::cerr << m_className << "::OptimiseOnTrack: "
10820 << "Error in the field function.\n";
10821 return false;
10822 }
10823 typedef std::function<double(const std::vector<double>&)> Fcn;
10824 Fcn& fField = *((Fcn *) fPtr);
10825
10826 // TODO: add option to use the current average of the field function
10827 // as target value.
10828
10829 // Convert the target function.
10830 if (m_polar) {
10831 variables = {"r", "phi"};
10832 } else {
10833 variables = {"x", "y"};
10834 }
10835 fPtr = MakeFunction("fTarget", target_function, variables);
10836 if (!fPtr) {
10837 std::cerr << m_className << "::OptimiseOnTrack: "
10838 << "Error in the target function.\n";
10839 return false;
10840 }
10841 Fcn& fTarget = *((Fcn *) fPtr);
10842
10843 // TODO: weight function.
10844
10845 std::vector<double> xFit(nP);
10846 std::iota(xFit.begin(), xFit.end(), 0);
10847 std::vector<double> yFit(nP);
10848 std::vector<double> wFit(nP, 1.);
10849 const double dx = (x1 - x0) / double(nP - 1);
10850 const double dy = (y1 - y0) / double(nP - 1);
10851 // OPTXYA
10852 for (unsigned int i = 0; i < nP; ++i) {
10853 // Position variables.
10854 std::vector<double> var = {x0 + i * dx, y0 + i * dy};
10855 if (m_polar) Cartesian2Polar(var[0], var[1], var[0], var[1]);
10856 // Evaluate the position dependent target function.
10857 yFit[i] = fTarget(var);
10858 // Evaluate the position dependent weighting function.
10859 // TODO.
10860 }
10861
10862 // Initialise the fit parameters.
10863 std::vector<double> vw0;
10864 std::array<double, 5> vp0;
10865 std::vector<double> aFit;
10866 std::vector<std::vector<unsigned int> > wiresInGroup;
10867 std::vector<std::vector<unsigned int> > planesInGroup;
10868 InitialiseFitParameters(groups, vw0, vp0, aFit, wiresInGroup, planesInGroup);
10869 if (aFit.empty()) {
10870 std::cerr << m_className << "::OptimiseOnTrack: "
10871 << "Setting fitting parameters failed.\n";
10872 return false;
10873 }
10874
10875 auto fFit = [&](const double s, const std::vector<double>& par) {
10876 // -------------------------------------------------------------------
10877 // OPTFUN - Function returning the value of the field function at a
10878 // given position (integer code).
10879 //--------------------------------------------------------------------
10880 // First set the potentials.
10881 std::vector<double> vw = vw0;
10882 std::array<double, 5> vp = vp0;
10883 const size_t nPar = par.size();
10884 for (size_t i = 0; i < nPar; ++i) {
10885 for (unsigned int iw : wiresInGroup[i]) vw[iw] += par[i];
10886 for (size_t ip : planesInGroup[i]) vp[ip] += par[i];
10887 }
10888 // Next reconstruct the charges.
10889 if (!Update(vw, vp)) return 0.;
10890
10891 std::vector<double> var(6, 0.);
10892 var[0] = x0 + s * dx;
10893 var[1] = y0 + s * dy;
10894 if (m_polar) Cartesian2Internal(var[0], var[1], var[0], var[1]);
10895 double ez = 0.;
10896 Field(var[0], var[1], 0., var[2], var[3], ez, var[5], true);
10897 var[4] = sqrt(var[2] * var[2] + var[3] * var[3]);
10898 // TODO: drift time, diffusion, gain, ...
10899 // Transform to polar if needed.
10900 if (m_polar) {
10901 Internal2Polar(var[0], var[1], var[0], var[1]);
10902 var[2] = var[2] / var[0];
10903 var[3] = var[3] / var[0];
10904 var[4] = var[4] / var[0];
10905 }
10906 // Calculate the field function with this field.
10907 return fField(var);
10908 };
10909
10910 double chi2 = 0.;
10911 std::vector<double> eFit(aFit.size(), 0.);
10912 // Carry out the fitting itself.
10913 if (!Numerics::LeastSquaresFit(fFit, aFit, eFit, xFit, yFit, wFit,
10914 m_optNitmax, m_optDist, chi2, m_optEps,
10915 print, m_debug)) {
10916 std::cerr << m_className << "::OptimiseOnTrack:\n"
10917 << " The new potentials do not fulfill your requirements.\n";
10918 return false;
10919 }
10920
10921 // Calculate the charges for the final result.
10922 std::vector<double> vw = vw0;
10923 std::array<double, 5> vp = vp0;
10924 for (size_t i = 0; i < aFit.size(); ++i) {
10925 for (unsigned int iw : wiresInGroup[i]) vw[iw] += aFit[i];
10926 for (size_t ip : planesInGroup[i]) vp[ip] += aFit[i];
10927 }
10928 if (!Update(vw, vp)) {
10929 std::cerr << m_className << "::OptimiseOnTrack:\n"
10930 << " Failed to compute the charges for the final settings.\n";
10931 return false;
10932 }
10933 return true;
10934}

◆ OptimiseOnWires()

bool Garfield::ComponentAnalyticField::OptimiseOnWires ( const std::vector< std::string > & groups,
const std::string & field_function,
const double target,
const std::vector< unsigned int > & wires,
const bool print = true )

Vary the potential of selected electrodes to match (a function of) the potential or field on the surfaces of a set of wires.

Definition at line 11093 of file ComponentAnalyticField.cc.

11096 {
11097 // -----------------------------------------------------------------------
11098 // OPTSET - Routine attempting to find proper voltage settings.
11099 // -----------------------------------------------------------------------
11100
11101 if (wires.empty()) {
11102 std::cerr << "OptimiseOnWires: Number of wires < 1.\n";
11103 return false;
11104 }
11105
11106 if (groups.empty()) {
11107 std::cerr << m_className << "::OptimiseOnWires: No electrode groups.\n";
11108 return false;
11109 }
11110
11111 // Use a constant target value for now.
11112 std::string target_function = std::to_string(target);
11113
11114 if (m_debug) {
11115 std::cout << m_className << "::OptimiseOnWires:\n"
11116 << " The function " << field_function << "\n"
11117 << " has to approximate the value of the function "
11118 << target_function << ".\n";
11119 }
11120
11121 // Convert the field function.
11122 std::vector<std::string> variables;
11123 if (m_polar) {
11124 variables = {"r", "phi", "er", "ephi", "e", "v"};
11125 } else {
11126 variables = {"x", "y", "ex", "ey", "e", "v"};
11127 }
11128 auto fPtr = MakeFunction("fField", field_function, variables);
11129 if (!fPtr) {
11130 std::cerr << m_className << "::OptimiseOnWires: "
11131 << "Error in the field function.\n";
11132 return false;
11133 }
11134 typedef std::function<double(const std::vector<double>&)> Fcn;
11135 Fcn& fField = *((Fcn *) fPtr);
11136
11137 // TODO: add option to use the current average of the field function
11138 // as target value.
11139
11140 // Convert the target function.
11141 if (m_polar) {
11142 variables = {"r", "phi"};
11143 } else {
11144 variables = {"x", "y"};
11145 }
11146 fPtr = MakeFunction("fTarget", target_function, variables);
11147 if (!fPtr) {
11148 std::cerr << m_className << "::OptimiseOnWires: "
11149 << "Error in the target function.\n";
11150 return false;
11151 }
11152 Fcn& fTarget = *((Fcn *) fPtr);
11153
11154 const size_t nW = wires.size();
11155 std::vector<double> xFit(nW, 0.);
11156 std::iota(xFit.begin(), xFit.end(), 0);
11157 std::vector<double> wFit(nW, 1.);
11158 std::vector<double> yFit;
11159 for (unsigned int iw : wires) {
11160 if (iw >= m_nWires) {
11161 std::cerr << m_className << "::OptimiseOnWires:\n"
11162 << " Wire index " << iw << " out of range.\n";
11163 return false;
11164 }
11165 // Position.
11166 std::vector<double> var = {m_w[iw].x, m_w[iw].y};
11167 if (m_polar) Internal2Polar(var[0], var[1], var[0], var[1]);
11168 // Evaluate the position dependent target function.
11169 yFit.push_back(fTarget(var));
11170 // Evaluate the position dependent weighting function.
11171 // TODO
11172 }
11173
11174 // Initialise the fit parameters.
11175 std::vector<double> vw0;
11176 std::array<double, 5> vp0;
11177 std::vector<double> aFit;
11178 std::vector<std::vector<unsigned int> > wiresInGroup;
11179 std::vector<std::vector<unsigned int> > planesInGroup;
11180 InitialiseFitParameters(groups, vw0, vp0, aFit, wiresInGroup, planesInGroup);
11181 if (aFit.empty()) {
11182 std::cerr << m_className << "::OptimiseOnWires: "
11183 << "Setting fitting parameters failed.\n";
11184 return false;
11185 }
11186
11187 auto fFit = [&](const double s, const std::vector<double>& par) {
11188 // First set the potentials.
11189 std::vector<double> vw = vw0;
11190 std::array<double, 5> vp = vp0;
11191 const size_t nPar = par.size();
11192 for (size_t i = 0; i < nPar; ++i) {
11193 for (unsigned int iw : wiresInGroup[i]) vw[iw] += par[i];
11194 for (size_t ip : planesInGroup[i]) vp[ip] += par[i];
11195 }
11196 // Next reconstruct the charges.
11197 if (!Update(vw, vp)) return 0.;
11198
11199 std::vector<double> var(6, 0.);
11200 const unsigned int iw = wires[int(s)];
11201 const double rw = m_w[iw].r;
11202 const double xw = m_w[iw].x;
11203 const double yw = m_w[iw].y;
11204 m_w[iw].r = 0.;
11205 var[0] = xw;
11206 var[1] = yw;
11207 constexpr unsigned int nA = 20;
11208 constexpr double dphi = TwoPi / nA;
11209 for (unsigned int i = 0; i < nA; ++i) {
11210 double phi = i * dphi;
11211 double ex = 0., ey = 0., ez = 0., volt = 0.;
11212 Field(xw + cos(phi) * rw, yw + sin(phi) * rw, 0.,
11213 ex, ey, ez, volt, true);
11214 var[4] += sqrt(ex * ex + ey * ey);
11215 var[5] += volt;
11216 }
11217 var[4] /= nA;
11218 var[5] /= nA;
11219 m_w[iw].r = rw;
11220 // Transform to polar if needed.
11221 if (m_polar) {
11222 Internal2Polar(var[0], var[1], var[0], var[1]);
11223 var[2] = var[2] / var[0];
11224 var[3] = var[3] / var[0];
11225 var[4] = var[4] / var[0];
11226 }
11227 // Calculate the field function with this field.
11228 return fField(var);
11229 };
11230
11231 double chi2 = 0.;
11232 std::vector<double> eFit(aFit.size(), 0.);
11233 // Carry out the fitting itself.
11234 if (!Numerics::LeastSquaresFit(fFit, aFit, eFit, xFit, yFit, wFit,
11235 m_optNitmax, m_optDist, chi2, m_optEps,
11236 print, m_debug)) {
11237 std::cerr << m_className << "::OptimiseOnWires:\n"
11238 << " The new potentials do not fulfill your requirements.\n";
11239 return false;
11240 }
11241
11242 // Calculate the charges for the final result.
11243 std::vector<double> vw = vw0;
11244 std::array<double, 5> vp = vp0;
11245 for (size_t i = 0; i < aFit.size(); ++i) {
11246 for (unsigned int iw : wiresInGroup[i]) vw[iw] += aFit[i];
11247 for (size_t ip : planesInGroup[i]) vp[ip] += aFit[i];
11248 }
11249 if (!Update(vw, vp)) {
11250 std::cerr << m_className << "::OptimiseOnWires:\n"
11251 << " Failed to compute the charges for the final settings.\n";
11252 return false;
11253 }
11254 return true;
11255}

◆ PlotCell()

void Garfield::ComponentAnalyticField::PlotCell ( TPad * pad)

Make a plot of the cell layout.

Definition at line 690 of file ComponentAnalyticField.cc.

690 {
691
692 ViewCell view;
693 view.SetComponent(this);
694 if (pad) view.SetCanvas(pad);
695 view.Plot2d();
696}

◆ PrintCell()

void Garfield::ComponentAnalyticField::PrintCell ( )

Print all available information on the cell.

Definition at line 387 of file ComponentAnalyticField.cc.

387 {
388 //-----------------------------------------------------------------------
389 // CELPRT - Subroutine printing all available information on the cell.
390 //-----------------------------------------------------------------------
391
392 // Make sure the cell is prepared.
393 if (!m_cellset && !Prepare()) {
394 std::cerr << m_className << "::PrintCell: Cell not set up.\n";
395 return;
396 }
397 std::cout << m_className
398 << "::PrintCell: Cell identification: " << GetCellType() << "\n";
399 // Print positions of wires, applied voltages and resulting charges.
400 if (!m_w.empty()) {
401 std::cout << " Table of the wires\n";
402 if (m_polar) {
403 std::cout << " Nr Diameter r phi Voltage";
404 } else {
405 std::cout << " Nr Diameter x y Voltage";
406 }
407 std::cout << " Charge Tension Length Density Label\n";
408 if (m_polar) {
409 std::cout << " [micron] [cm] [deg] [Volt]";
410 } else {
411 std::cout << " [micron] [cm] [cm] [Volt]";
412 }
413 std::cout << " [pC/cm] [g] [cm] [g/cm3]\n";
414 for (unsigned int i = 0; i < m_nWires; ++i) {
415 const auto& w = m_w[i];
416 double xw = w.x;
417 double yw = w.y;
418 double dw = 2 * w.r;
419 if (m_polar) {
420 Internal2Polar(w.x, w.y, xw, yw);
421 dw *= xw;
422 }
423 std::printf(
424 "%4u %9.2f %9.4f %9.4f %9.3f %12.4f %9.2f %9.2f %9.2f \"%s\"\n", i,
425 1.e4 * dw, xw, yw, w.v, w.e * TwoPiEpsilon0 * 1.e-3, w.tension,
426 w.u, w.density, w.type.c_str());
427 }
428 }
429 // Print information on the tube if present.
430 if (m_tube) {
431 std::string shape;
432 if (m_ntube == 0) {
433 shape = "Circular";
434 } else if (m_ntube == 3) {
435 shape = "Triangular";
436 } else if (m_ntube == 4) {
437 shape = "Square";
438 } else if (m_ntube == 5) {
439 shape = "Pentagonal";
440 } else if (m_ntube == 6) {
441 shape = "Hexagonal";
442 } else if (m_ntube == 7) {
443 shape = "Heptagonal";
444 } else if (m_ntube == 8) {
445 shape = "Octagonal";
446 } else {
447 shape = "Polygonal with " + std::to_string(m_ntube) + " corners";
448 }
449 std::cout << " Enclosing tube\n"
450 << " Potential: " << m_vttube << " V\n"
451 << " Radius: " << m_cotube << " cm\n"
452 << " Shape: " << shape << "\n"
453 << " Label: " << m_planes[4].type << "\n";
454 }
455 // Print data on the equipotential planes.
456 if (m_ynplan[0] || m_ynplan[1] || m_ynplan[2] || m_ynplan[3]) {
457 std::cout << " Equipotential planes\n";
458 // First those at const x or r.
459 const std::string xr = m_polar ? "r" : "x";
460 if (m_ynplan[0] && m_ynplan[1]) {
461 std::cout << " There are two planes at constant " << xr << ":\n";
462 } else if (m_ynplan[0] || m_ynplan[1]) {
463 std::cout << " There is one plane at constant " << xr << ":\n";
464 }
465 for (unsigned int i = 0; i < 2; ++i) {
466 if (!m_ynplan[i]) continue;
467 if (m_polar) {
468 std::cout << " r = " << exp(m_coplan[i]) << " cm, ";
469 } else {
470 std::cout << " x = " << m_coplan[i] << " cm, ";
471 }
472 if (fabs(m_vtplan[i]) > 1.e-4) {
473 std::cout << "potential = " << m_vtplan[i] << " V, ";
474 } else {
475 std::cout << "earthed, ";
476 }
477 const auto& plane = m_planes[i];
478 if (!plane.type.empty() && plane.type != "?") {
479 std::cout << "label = " << plane.type << ", ";
480 }
481 const unsigned int nStrips = plane.strips1.size() + plane.strips2.size();
482 const unsigned int nPixels = plane.pixels.size();
483 if (nStrips == 0 && nPixels == 0) {
484 std::cout << "no strips or pixels.\n";
485 } else if (nPixels == 0) {
486 if (nStrips == 1) std::cout << "1 strip.\n";
487 else std::cout << nStrips << " strips.\n";
488 } else if (nStrips == 0) {
489 if (nPixels == 1) std::cout << "1 pixel.\n";
490 else std::cout << nPixels << " pixels.\n";
491 } else {
492 if (nStrips == 1) std::cout << "1 strip, ";
493 else std::cout << nStrips << " strips, ";
494 if (nPixels == 1) std::cout << "1 pixel.\n";
495 else std::cout << nPixels << " pixels.\n";
496 }
497 for (const auto& strip : plane.strips2) {
498 std::cout << " ";
499 if (m_polar) {
500 double gap = i == 0 ? expm1(strip.gap) : -expm1(-strip.gap);
501 gap *= exp(m_coplan[i]);
502 std::cout << RadToDegree * strip.smin << " < phi < "
503 << RadToDegree * strip.smax
504 << " degrees, gap = " << gap << " cm";
505 } else {
506 std::cout << strip.smin << " < y < " << strip.smax
507 << " cm, gap = " << strip.gap << " cm";
508 }
509 if (!strip.type.empty() && strip.type != "?") {
510 std::cout << " (\"" << strip.type << "\")";
511 }
512 std::cout << "\n";
513 }
514 for (const auto& strip : plane.strips1) {
515 std::cout << " " << strip.smin << " < z < " << strip.smax;
516 if (m_polar) {
517 double gap = i == 0 ? expm1(strip.gap) : -expm1(-strip.gap);
518 gap *= exp(m_coplan[i]);
519 std::cout << " cm, gap = " << gap << " cm";
520 } else {
521 std::cout << " cm, gap = " << strip.gap << " cm";
522 }
523 if (!strip.type.empty() && strip.type != "?") {
524 std::cout << " (\"" << strip.type << "\")";
525 }
526 std::cout << "\n";
527 }
528 for (const auto& pix : plane.pixels) {
529 std::cout << " ";
530 if (m_polar) {
531 std::cout << RadToDegree * pix.smin << " < phi < "
532 << RadToDegree * pix.smax << " degrees, ";
533 } else {
534 std::cout << pix.smin << " < y < " << pix.smax << " cm, ";
535 }
536 std::cout << pix.zmin << " < z < " << pix.zmax << " cm, gap = ";
537 if (m_polar) {
538 double gap = i == 0 ? expm1(pix.gap) : -expm1(-pix.gap);
539 gap *= exp(m_coplan[i]);
540 std::cout << gap << " cm";
541 } else {
542 std::cout << pix.gap << " cm";
543 }
544 if (!pix.type.empty() && pix.type != "?") {
545 std::cout << " (\"" << pix.type << "\")";
546 }
547 std::cout << "\n";
548 }
549 }
550 // Next the planes at constant y or phi
551 const std::string yphi = m_polar ? "phi" : "y";
552 if (m_ynplan[2] && m_ynplan[3]) {
553 std::cout << " There are two planes at constant " << yphi << ":\n";
554 } else if (m_ynplan[2] || m_ynplan[3]) {
555 std::cout << " There is one plane at constant " << yphi << ":\n";
556 }
557 for (unsigned int i = 2; i < 4; ++i) {
558 if (!m_ynplan[i]) continue;
559 if (m_polar) {
560 std::cout << " phi = " << RadToDegree * m_coplan[i] << " degrees, ";
561 } else {
562 std::cout << " y = " << m_coplan[i] << " cm, ";
563 }
564 if (fabs(m_vtplan[i]) > 1.e-4) {
565 std::cout << "potential = " << m_vtplan[i] << " V, ";
566 } else {
567 std::cout << "earthed, ";
568 }
569 const auto& plane = m_planes[i];
570 if (!plane.type.empty() && plane.type != "?") {
571 std::cout << "label = " << plane.type << ", ";
572 }
573 const unsigned int nStrips = plane.strips1.size() + plane.strips2.size();
574 const unsigned int nPixels = plane.pixels.size();
575 if (nStrips == 0 && nPixels == 0) {
576 std::cout << "no strips or pixels.\n";
577 } else if (nPixels == 0) {
578 std::cout << nStrips << " strips.\n";
579 } else if (nStrips == 0) {
580 std::cout << nPixels << " pixels.\n";
581 } else {
582 std::cout << nStrips << " strips, " << nPixels << " pixels.\n";
583 }
584 for (const auto& strip : plane.strips2) {
585 std::cout << " ";
586 if (m_polar) {
587 std::cout << exp(strip.smin) << " < r < " << exp(strip.smax)
588 << " cm, gap = " << RadToDegree * strip.gap << " degrees";
589 } else {
590 std::cout << strip.smin << " < x < " << strip.smax
591 << " cm, gap = " << strip.gap << " cm";
592 }
593 if (!strip.type.empty() && strip.type != "?") {
594 std::cout << " (\"" << strip.type << "\")";
595 }
596 std::cout << "\n";
597 }
598 for (const auto& strip : plane.strips1) {
599 std::cout << " " << strip.smin << " < z < " << strip.smax;
600 if (m_polar) {
601 std::cout << " cm, gap = " << RadToDegree * strip.gap << " degrees";
602 } else {
603 std::cout << " cm, gap = " << strip.gap << " cm";
604 }
605 if (!strip.type.empty() && strip.type != "?") {
606 std::cout << " (\"" << strip.type << "\")";
607 }
608 std::cout << "\n";
609 }
610 for (const auto& pix : plane.pixels) {
611 std::cout << " ";
612 if (m_polar) {
613 std::cout << exp(pix.smin) << " < r < " << exp(pix.smax) << " cm, ";
614 } else {
615 std::cout << pix.smin << " < x < " << pix.smax << " cm, ";
616 }
617 std::cout << pix.zmin << " < z < " << pix.zmax << " cm, gap = ";
618 if (m_polar) {
619 std::cout << RadToDegree * pix.gap << " degrees";
620 } else {
621 std::cout << pix.gap << " cm";
622 }
623 if (!pix.type.empty() && pix.type != "?") {
624 std::cout << " (\"" << pix.type << "\")";
625 }
626 std::cout << "\n";
627 }
628 }
629 }
630 // Print the type of periodicity.
631 std::cout << " Periodicity\n";
632 if (m_perx) {
633 std::cout << " The cell is repeated every ";
634 if (m_polar) {
635 std::cout << exp(m_sx) << " cm in r.\n";
636 } else {
637 std::cout << m_sx << " cm in x.\n";
638 }
639 } else {
640 if (m_polar) {
641 std::cout << " The cell is not periodic in r.\n";
642 } else {
643 std::cout << " The cell has no translation periodicity in x.\n";
644 }
645 }
646 if (m_pery) {
647 std::cout << " The cell is repeated every ";
648 if (m_polar) {
649 std::cout << RadToDegree * m_sy << " degrees in phi.\n";
650 } else {
651 std::cout << m_sy << " cm in y.\n";
652 }
653 } else {
654 if (m_polar) {
655 std::cout << " The cell is not periodic in phi.\n";
656 } else {
657 std::cout << " The cell has no translation periodicity in y.\n";
658 }
659 }
660 std::cout << " Other data\n";
661 std::cout << " Gravity vector: (" << m_down[0] << ", " << m_down[1]
662 << ", " << m_down[2] << ") [g].\n";
663 std::cout << " Cell dimensions:\n ";
664 if (!m_polar) {
665 std::cout << m_xmin << " < x < " << m_xmax << " cm, " << m_ymin << " < y < "
666 << m_ymax << " cm.\n";
667 } else {
668 double xminp, yminp;
669 Internal2Polar(m_xmin, m_ymin, xminp, yminp);
670 double xmaxp, ymaxp;
671 Internal2Polar(m_xmax, m_ymax, xmaxp, ymaxp);
672 std::cout << xminp << " < r < " << xmaxp << " cm, " << yminp << " < phi < "
673 << ymaxp << " degrees.\n";
674 }
675 std::cout << " Potential range:\n " << m_vmin << " < V < " << m_vmax
676 << " V.\n";
677 // Print voltage shift in case no equipotential planes are present.
678 if (!(m_ynplan[0] || m_ynplan[1] || m_ynplan[2] || m_ynplan[3] || m_tube)) {
679 std::cout << " All voltages have been shifted by " << m_v0
680 << " V to avoid net wire charge.\n";
681 } else {
682 // Print the net charge on wires.
683 double sum = 0.;
684 for (const auto& w : m_w) sum += w.e;
685 std::cout << " The net charge on the wires is "
686 << 1.e-3 * TwoPiEpsilon0 * sum << " pC/cm.\n";
687 }
688}

◆ PrintCharges()

void Garfield::ComponentAnalyticField::PrintCharges ( ) const

Print a list of the point charges.

Definition at line 1781 of file ComponentAnalyticField.cc.

1781 {
1782 std::cout << m_className << "::PrintCharges:\n";
1783 if (m_ch3d.empty()) {
1784 std::cout << " No charges present.\n";
1785 return;
1786 }
1787 std::cout << " x [cm] y [cm] z [cm] charge [fC]\n";
1788 for (const auto& charge : m_ch3d) {
1789 std::cout << " " << std::setw(9) << charge.x << " " << std::setw(9)
1790 << charge.y << " " << std::setw(9) << charge.z << " "
1791 << std::setw(11) << charge.e * FourPiEpsilon0 << "\n";
1792 }
1793}

◆ SetCartesianCoordinates()

void Garfield::ComponentAnalyticField::SetCartesianCoordinates ( )

Use Cartesian coordinates (default).

Definition at line 1743 of file ComponentAnalyticField.cc.

1743 {
1744 if (m_polar) {
1745 std::cout << m_className << "::SetCartesianCoordinates:\n "
1746 << "Switching to Cartesian coordinates; resetting the cell.\n";
1747 CellInit();
1748 }
1749 m_polar = false;
1750}

◆ SetGravity()

void Garfield::ComponentAnalyticField::SetGravity ( const double dx,
const double dy,
const double dz )

Set the gravity orientation.

Definition at line 2055 of file ComponentAnalyticField.cc.

2056 {
2057
2058 const double d = sqrt(dx * dx + dy * dy + dz * dz);
2059 if (d > 0.) {
2060 m_down[0] = dx / d;
2061 m_down[1] = dy / d;
2062 m_down[2] = dz / d;
2063 } else {
2064 std::cerr << m_className << "::SetGravity:\n"
2065 << " The gravity vector has zero norm ; ignored.\n";
2066 }
2067}

◆ SetMedium()

void Garfield::ComponentAnalyticField::SetMedium ( Medium * medium)
inline

Set the medium inside the cell.

Definition at line 26 of file ComponentAnalyticField.hh.

26{ m_medium = medium; }

Referenced by main().

◆ SetNumberOfCellCopies()

void Garfield::ComponentAnalyticField::SetNumberOfCellCopies ( const unsigned int nfourier)

Definition at line 3837 of file ComponentAnalyticField.cc.

3837 {
3838 if (nf > 0) {
3839 if ((nf & (nf - 1)) != 0) {
3840 std::cerr << m_className << "::SetNumberOfCellCopies:\n"
3841 << " Argument must be an integral power of 2.\n";
3842 return;
3843 }
3844 }
3845 m_nFourier = nf;
3846 m_sigset = false;
3847}

◆ SetNumberOfShots()

void Garfield::ComponentAnalyticField::SetNumberOfShots ( const unsigned int n)
inline

Set the number of shots used for numerically solving the wire sag differential equation.

Definition at line 320 of file ComponentAnalyticField.hh.

320{ m_nShots = n; }

◆ SetNumberOfSteps()

void Garfield::ComponentAnalyticField::SetNumberOfSteps ( const unsigned int n)

Set the number of integration steps within each shot (must be >= 1).

Definition at line 2286 of file ComponentAnalyticField.cc.

2286 {
2287 if (n == 0) {
2288 std::cerr << m_className << "::SetNumberOfSteps:\n"
2289 << " Number of steps must be > 0.\n";
2290 return;
2291 }
2292 m_nSteps = n;
2293}

◆ SetOptimisationParameters()

void Garfield::ComponentAnalyticField::SetOptimisationParameters ( const double dist = 1.,
const double eps = 1.e-4,
const unsigned int nMaxIter = 10 )

Set the conditions at which to allow the iteration to stop.

Parameters
distMaximum deviation among all points between target and field function.
epsRelative change between iterations.
nMaxIterMaximum number of iterations.

Definition at line 10757 of file ComponentAnalyticField.cc.

10758 {
10759
10760 if (dist <= 0.) {
10761 std::cerr << m_className << "::SetOptimisationParameters: "
10762 << "Threshold distance must be > 0.\n";
10763 } else {
10764 m_optDist = dist;
10765 }
10766 if (eps <= 0.) {
10767 std::cerr << m_className << "::SetOptimisationParameters: "
10768 << "Relative change (epsilon) must be >= 0.\n";
10769 } else {
10770 m_optEps = eps;
10771 }
10772 if (iterlim < 1) {
10773 std::cerr << m_className << "::SetOptimisationParameters: "
10774 << "Maximum number of iterations must be > 0.\n";
10775 } else {
10776 m_optNitmax = iterlim;
10777 }
10778}

◆ SetPeriodicityPhi()

void Garfield::ComponentAnalyticField::SetPeriodicityPhi ( const double phi)

Set the periodicity [degree] in phi.

Definition at line 1638 of file ComponentAnalyticField.cc.

1638 {
1639 if (!m_polar && !m_tube) {
1640 std::cerr << m_className << "::SetPeriodicityPhi:\n"
1641 << " Cannot use phi-periodicity with Cartesian coordinates.\n";
1642 return;
1643 }
1644 if (std::abs(360. - s * int(360. / s)) > 1.e-4) {
1645 std::cerr << m_className << "::SetPeriodicityPhi:\n"
1646 << " Phi periods must divide 360; ignored.\n";
1647 return;
1648 }
1649
1650 m_periodic[1] = true;
1651 m_sy = DegreeToRad * s;
1652 m_mtube = int(360. / s);
1653 UpdatePeriodicity();
1654}

◆ SetPeriodicityX()

void Garfield::ComponentAnalyticField::SetPeriodicityX ( const double s)

Set the periodic length [cm] in the x-direction.

Definition at line 1604 of file ComponentAnalyticField.cc.

1604 {
1605 if (m_polar) {
1606 std::cerr << m_className << "::SetPeriodicityX:\n"
1607 << " Cannot use x-periodicity with polar coordinates.\n";
1608 return;
1609 }
1610 if (s < Small) {
1611 std::cerr << m_className << "::SetPeriodicityX:\n"
1612 << " Periodic length must be greater than zero.\n";
1613 return;
1614 }
1615
1616 m_periodic[0] = true;
1617 m_sx = s;
1618 UpdatePeriodicity();
1619}

◆ SetPeriodicityY()

void Garfield::ComponentAnalyticField::SetPeriodicityY ( const double s)

Set the periodic length [cm] in the y-direction.

Definition at line 1621 of file ComponentAnalyticField.cc.

1621 {
1622 if (m_polar) {
1623 std::cerr << m_className << "::SetPeriodicityY:\n"
1624 << " Cannot use y-periodicity with polar coordinates.\n";
1625 return;
1626 }
1627 if (s < Small) {
1628 std::cerr << m_className << "::SetPeriodicityY:\n"
1629 << " Periodic length must be greater than zero.\n";
1630 return;
1631 }
1632
1633 m_periodic[1] = true;
1634 m_sy = s;
1635 UpdatePeriodicity();
1636}

◆ SetPolarCoordinates()

void Garfield::ComponentAnalyticField::SetPolarCoordinates ( )

Use polar coordinates.

Definition at line 1752 of file ComponentAnalyticField.cc.

1752 {
1753 if (!m_polar) {
1754 std::cout << m_className << "::SetPolarCoordinates:\n "
1755 << "Switching to polar coordinates; resetting the cell.\n";
1756 CellInit();
1757 }
1758 m_polar = true;
1759 // Set default phi period.
1760 m_pery = true;
1761 m_sy = TwoPi;
1762}

◆ SetScanningArea()

void Garfield::ComponentAnalyticField::SetScanningArea ( const double xmin,
const double xmax,
const double ymin,
const double ymax )

Specify explicitly the boundaries of the the scanning area (i. e. the area in which the electrostatic force acting on a wire is computed).

Definition at line 2029 of file ComponentAnalyticField.cc.

2032 {
2033 if (std::abs(xmax - xmin) < Small || std::abs(ymax - ymin) < Small) {
2034 std::cerr << m_className << "::SetScanningArea:\n"
2035 << " Zero range not permitted.\n";
2036 return;
2037 }
2038 m_scanRange = ScanningRange::User;
2039 m_xScanMin = std::min(xmin, xmax);
2040 m_xScanMax = std::max(xmin, xmax);
2041 m_yScanMin = std::min(ymin, ymax);
2042 m_yScanMax = std::max(ymin, ymax);
2043}

◆ SetScanningAreaFirstOrder()

void Garfield::ComponentAnalyticField::SetScanningAreaFirstOrder ( const double scale = 2.)

Set the scanning area based on the zeroth-order estimates of the wire shift, enlarged by a scaling factor. This is the default behaviour.

Definition at line 2045 of file ComponentAnalyticField.cc.

2045 {
2046 m_scanRange = ScanningRange::FirstOrder;
2047 if (scale > 0.) {
2048 m_scaleRange = scale;
2049 } else {
2050 std::cerr << m_className << "::SetScanningAreaFirstOrder:\n"
2051 << " Scaling factor must be > 0.\n";
2052 }
2053}

◆ SetScanningAreaLargest()

void Garfield::ComponentAnalyticField::SetScanningAreaLargest ( )
inline

Set the scanning area to the largest area around each wire which is free from other cell elements.

Definition at line 275 of file ComponentAnalyticField.hh.

275{ m_scanRange = ScanningRange::Largest; }

◆ SetScanningGrid()

void Garfield::ComponentAnalyticField::SetScanningGrid ( const unsigned int nX,
const unsigned int nY )

Set the number of grid lines at which the electrostatic force as function of the wire displacement is computed.

Definition at line 2013 of file ComponentAnalyticField.cc.

2014 {
2015 if (nX < 2) {
2016 std::cerr << m_className << "::SetScanningGrid:\n"
2017 << " Number of x-lines must be > 1.\n";
2018 } else {
2019 m_nScanX = nX;
2020 }
2021 if (nY < 2) {
2022 std::cerr << m_className << "::SetScanningGrid:\n"
2023 << " Number of y-lines must be > 1.\n";
2024 } else {
2025 m_nScanY = nY;
2026 }
2027}

◆ StepSizeHint()

double Garfield::ComponentAnalyticField::StepSizeHint ( )
overridevirtual

Reimplemented from Garfield::Component.

Definition at line 381 of file ComponentAnalyticField.cc.

381 {
382
383 if (!m_cellset && !Prepare()) return -1.;
384 return m_dmin;
385}

◆ WeightingField()

void Garfield::ComponentAnalyticField::WeightingField ( const double x,
const double y,
const double z,
double & wx,
double & wy,
double & wz,
const std::string & label )
inlineoverridevirtual

Calculate the weighting field at a given point and for a given electrode.

Parameters
x,y,zcoordinates [cm].
wx,wy,wzcomponents of the weighting field [1/cm].
labelname of the electrode

Reimplemented from Garfield::Component.

Definition at line 361 of file ComponentAnalyticField.hh.

363 {
364 wx = wy = wz = 0.;
365 if (!m_sigset) PrepareSignals();
366 Wfield(x, y, z, wx, wy, wz, label);
367 }

◆ WeightingPotential()

double Garfield::ComponentAnalyticField::WeightingPotential ( const double x,
const double y,
const double z,
const std::string & label )
inlineoverridevirtual

Calculate the weighting potential at a given point.

Parameters
x,y,zcoordinates [cm].
labelname of the electrode.
Returns
weighting potential [dimensionless].

Reimplemented from Garfield::Component.

Definition at line 368 of file ComponentAnalyticField.hh.

369 {
370 if (!m_sigset) PrepareSignals();
371 return Wpot(x, y, z, label);
372 }

◆ WireDisplacement()

bool Garfield::ComponentAnalyticField::WireDisplacement ( const unsigned int iw,
const bool detailed,
std::vector< double > & csag,
std::vector< double > & xsag,
std::vector< double > & ysag,
double & stretch,
const bool print = true )

Compute the sag profile of a wire.

Parameters
iwindex of the wire
detailedflag to request a detailed calculation of the profile or only a fast one.
csagcoordinate along the wire.
xsagx components of the sag profile.
ysagy components of the sag profile.
stretchrelative elongation.
printflag to print the calculation results or not.

Definition at line 2295 of file ComponentAnalyticField.cc.

2298 {
2299 if (!m_cellset && !Prepare()) {
2300 std::cerr << m_className << "::WireDisplacement: Cell not set up.\n";
2301 return false;
2302 }
2303 if (iw >= m_nWires) {
2304 std::cerr << m_className
2305 << "::WireDisplacement: Wire index out of range.\n";
2306 return false;
2307 }
2308 if (m_polar) {
2309 std::cerr << m_className
2310 << "::WireDisplacement: Cannot handle polar cells.\n";
2311 return false;
2312 }
2313 const auto& wire = m_w[iw];
2314 // Save the original coordinates.
2315 const double x0 = wire.x;
2316 const double y0 = wire.y;
2317 // First-order approximation.
2318 if (!Setup()) {
2319 std::cerr << m_className << "::WireDisplacement:\n"
2320 << " Charge calculation failed at central position.\n";
2321 return false;
2322 }
2323
2324 double fx = 0., fy = 0.;
2325 if (m_useElectrostaticForce) {
2326 double ex = 0., ey = 0.;
2327 ElectricFieldAtWire(iw, ex, ey);
2328 fx -= ex * wire.e * Internal2Newton;
2329 fy -= ey * wire.e * Internal2Newton;
2330 }
2331 if (m_useGravitationalForce) {
2332 // Mass per unit length [kg / cm].
2333 const double m = 1.e-3 * wire.density * Pi * wire.r * wire.r;
2334 fx -= m_down[0] * GravitationalAcceleration * m;
2335 fy -= m_down[1] * GravitationalAcceleration * m;
2336 }
2337 const double u2 = wire.u * wire.u;
2338 const double shiftx =
2339 -125 * fx * u2 / (GravitationalAcceleration * wire.tension);
2340 const double shifty =
2341 -125 * fy * u2 / (GravitationalAcceleration * wire.tension);
2342 // Get the elongation from this.
2343 const double s = 4 * sqrt(shiftx * shiftx + shifty * shifty) / wire.u;
2344 double length = wire.u;
2345 if (s > Small) {
2346 const double t = sqrt(1 + s * s);
2347 length *= 0.5 * (t + log(s + t) / s);
2348 }
2349 stretch = (length - wire.u) / wire.u;
2350 if (print) {
2351 std::cout << m_className
2352 << "::WireDisplacement:\n"
2353 << " Forces and displacement in 0th order.\n"
2354 << " Wire information: number = " << iw << "\n"
2355 << " type = " << wire.type << "\n"
2356 << " location = (" << wire.x << ", " << wire.y
2357 << ") cm\n"
2358 << " voltage = " << wire.v << " V\n"
2359 << " length = " << wire.u << " cm\n"
2360 << " tension = " << wire.tension << " g\n"
2361 << " In this position: Fx = " << fx << " N/cm\n"
2362 << " Fy = " << fy << " N/cm\n"
2363 << " x-shift = " << shiftx << " cm\n"
2364 << " y-shift = " << shifty << " cm\n"
2365 << " stretch = " << 100. * stretch << "%\n";
2366 }
2367 if (!detailed) {
2368 csag = {0.};
2369 xsag = {shiftx};
2370 ysag = {shifty};
2371 return true;
2372 }
2373 std::vector<double> xMap(m_nScanX, 0.);
2374 std::vector<double> yMap(m_nScanY, 0.);
2375 std::vector<std::vector<double> > fxMap(m_nScanX,
2376 std::vector<double>(m_nScanY, 0.));
2377 std::vector<std::vector<double> > fyMap(m_nScanX,
2378 std::vector<double>(m_nScanY, 0.));
2379 if (!ForcesOnWire(iw, xMap, yMap, fxMap, fyMap)) {
2380 std::cerr << m_className << "::WireDisplacement:\n"
2381 << " Could not compute the electrostatic force table.\n";
2382 return false;
2383 }
2384 // Compute the detailed wire shift.
2385 if (!SagDetailed(wire, xMap, yMap, fxMap, fyMap, csag, xsag, ysag)) {
2386 std::cerr << m_className << "::WireDisplacement: Wire " << iw << ".\n"
2387 << " Computation of the wire sag failed.\n";
2388 return false;
2389 }
2390 // Verify that the wire is in range.
2391 const double sxmin = xMap.front();
2392 const double sxmax = xMap.back();
2393 const double symin = yMap.front();
2394 const double symax = yMap.back();
2395 const unsigned int nSag = xsag.size();
2396 bool outside = false;
2397 length = 0.;
2398 double xAvg = 0.;
2399 double yAvg = 0.;
2400 double xMax = 0.;
2401 double yMax = 0.;
2402 for (unsigned int i = 0; i < nSag; ++i) {
2403 if (x0 + xsag[i] < sxmin || x0 + xsag[i] > sxmax || y0 + ysag[i] < symin ||
2404 y0 + ysag[i] > symax) {
2405 outside = true;
2406 }
2407 xAvg += xsag[i];
2408 yAvg += ysag[i];
2409 xMax = std::max(xMax, std::abs(xsag[i]));
2410 yMax = std::max(yMax, std::abs(ysag[i]));
2411 if (i == 0) continue;
2412 const double dx = xsag[i] - xsag[i - 1];
2413 const double dy = ysag[i] - ysag[i - 1];
2414 const double dz = csag[i] - csag[i - 1];
2415 length += sqrt(dx * dx + dy * dy + dz * dz);
2416 }
2417 xAvg /= nSag;
2418 yAvg /= nSag;
2419 stretch = (length - wire.u) / wire.u;
2420 // Warn if a point outside the scanning area was found.
2421 if (outside) {
2422 std::cerr
2423 << m_className << "::WireDisplacement: Warning.\n "
2424 << "The wire profile is located partially outside the scanning area.\n";
2425 }
2426 if (print) {
2427 std::cout << " Sag profile for wire " << iw << ".\n"
2428 << " Point z [cm] x-sag [um] y-sag [um]\n";
2429 for (unsigned int i = 0; i < nSag; ++i) {
2430 std::printf(" %3u %10.4f %10.4f %10.4f\n",
2431 i, csag[i], xsag[i] * 1.e4, ysag[i] * 1.e4);
2432 }
2433 std::printf(" Average sag in x and y: %10.4f and %10.4f micron\n",
2434 1.e4 * xAvg, 1.e4 * yAvg);
2435 std::printf(" Maximum sag in x and y: %10.4f and %10.4f micron\n",
2436 1.e4 * xMax, 1.e4 * yMax);
2437 std::cout << " Elongation: " << 100. * stretch << "%\n";
2438 }
2439 return true;
2440}
Garfield::ComponentAnalyticField::ForcesOnWire
bool ForcesOnWire(const unsigned int iw, std::vector< double > &xMap, std::vector< double > &yMap, std::vector< std::vector< double > > &fxMap, std::vector< std::vector< double > > &fyMap)
Definition ComponentAnalyticField.cc:2077
The documentation for this class was generated from the following files: