Garfield++ 3.0
A toolkit for the detailed simulation of particle detectors based on ionisation measurement in gases and semiconductors
Loading...
Searching...
No Matches
Classes | 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.
 
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 IsWireCrossed (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 IsInTrapRadius (const double q0, const double x0, const double y0, const double z0, double &xw, double &yx, double &rw) override
 
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.)
 Add a strip in the y or z direction on an existing plane at constant x.
 
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.)
 Add a pixel on an existing plane at constant x.
 
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.)
 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 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)
 Setup the weighting field for a given group of wires or planes.
 
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 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.
 
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).
 
- Public Member Functions inherited from Garfield::ComponentBase
 ComponentBase ()
 Constructor.
 
virtual ~ComponentBase ()
 Destructor.
 
virtual void SetGeometry (GeometryBase *geo)
 Define the geometry.
 
virtual void Clear ()
 Reset.
 
virtual Medium * GetMedium (const double x, const double y, const double z)
 Get the medium at a given location (x, y, z).
 
virtual void ElectricField (const double x, const double y, const double z, double &ex, double &ey, double &ez, Medium *&m, int &status)=0
 
virtual void ElectricField (const double x, const double y, const double z, double &ex, double &ey, double &ez, double &v, Medium *&m, int &status)=0
 Calculate the drift field [V/cm] and potential [V] at (x, y, z).
 
virtual bool GetVoltageRange (double &vmin, double &vmax)=0
 Calculate the voltage range [V].
 
virtual void WeightingField (const double x, const double y, const double z, double &wx, double &wy, double &wz, const std::string &label)
 
virtual double WeightingPotential (const double x, const double y, const double z, const std::string &label)
 
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 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?
 
virtual bool GetBoundingBox (double &xmin, double &ymin, double &zmin, double &xmax, double &ymax, double &zmax)
 Get the bounding box coordinates.
 
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 IntegrateFlux (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)
 
virtual bool IsWireCrossed (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)
 
virtual bool IsInTrapRadius (const double q0, const double x0, const double y0, const double z0, double &xw, double &yw, double &rw)
 
void EnablePeriodicityX (const bool on=true)
 Enable simple periodicity in the $x$ direction.
 
void DisablePeriodicityX ()
 
void EnablePeriodicityY (const bool on=true)
 Enable simple periodicity in the $y$ direction.
 
void DisablePeriodicityY ()
 
void EnablePeriodicityZ (const bool on=true)
 Enable simple periodicity in the $z$ direction.
 
void DisablePeriodicityZ ()
 
void EnableMirrorPeriodicityX (const bool on=true)
 Enable mirror periodicity in the $x$ direction.
 
void DisableMirrorPeriodicityX ()
 
void EnableMirrorPeriodicityY (const bool on=true)
 Enable mirror periodicity in the $y$ direction.
 
void DisableMirrorPeriodicityY ()
 
void EnableMirrorPeriodicityZ (const bool on=true)
 Enable mirror periodicity in the $y$ direction.
 
void DisableMirrorPeriodicityZ ()
 
void EnableAxialPeriodicityX (const bool on=true)
 Enable axial periodicity in the $x$ direction.
 
void DisableAxialPeriodicityX ()
 
void EnableAxialPeriodicityY (const bool on=true)
 Enable axial periodicity in the $y$ direction.
 
void DisableAxialPeriodicityY ()
 
void EnableAxialPeriodicityZ (const bool on=true)
 Enable axial periodicity in the $z$ direction.
 
void DisableAxialPeriodicityZ ()
 
void EnableRotationSymmetryX (const bool on=true)
 Enable rotation symmetry around the $x$ axis.
 
void DisableRotationSymmetryX ()
 
void EnableRotationSymmetryY (const bool on=true)
 Enable rotation symmetry around the $y$ axis.
 
void DisableRotationSymmetryY ()
 
void EnableRotationSymmetryZ (const bool on=true)
 Enable rotation symmetry around the $z$ axis.
 
void DisableRotationSymmetryZ ()
 
void EnableDebugging ()
 Switch on debugging messages.
 
void DisableDebugging ()
 Switch off debugging messages.
 
void ActivateTraps ()
 Request trapping to be taken care of by the component (for TCAD).
 
void DeactivateTraps ()
 
bool IsTrapActive ()
 
void ActivateVelocityMap ()
 Request velocity to be taken care of by the component (for TCAD).
 
void DectivateVelocityMap ()
 
bool IsVelocityActive ()
 
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 void ElectronVelocity (const double, const double, const double, double &vx, double &vy, double &vz, Medium *&, int &status)
 Get the electron drift velocity.
 
virtual void HoleVelocity (const double, const double, const double, double &vx, double &vy, double &vz, Medium *&, int &status)
 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

virtual void Reset ()=0
 Reset the component.
 
virtual void UpdatePeriodicity ()=0
 Verify periodicities.
 
- Protected Attributes inherited from Garfield::ComponentBase
std::string m_className = "ComponentBase"
 Class name.
 
GeometryBase * m_geometry = nullptr
 Pointer to the geometry.
 
bool m_ready = false
 Ready for use?
 
bool m_activeTraps = false
 Does the component have traps?
 
bool m_hasVelocityMap = false
 Does the component have velocity maps?
 
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.
 
double m_bx0 = 0.
 
double m_by0 = 0.
 
double m_bz0 = 0.
 
bool m_debug = false
 Switch on/off debugging messages.
 

Detailed Description

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

Definition at line 15 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 313 of file ComponentAnalyticField.hh.

313 {
314 A00,
315 B1X,
316 B1Y,
317 B2X,
318 B2Y,
319 C10,
320 C2X,
321 C2Y,
322 C30,
323 D10,
324 D20,
325 D30,
326 D40,
327 Unknown
328 };

Constructor & Destructor Documentation

◆ ComponentAnalyticField()

Garfield::ComponentAnalyticField::ComponentAnalyticField ( )

Constructor.

Definition at line 240 of file ComponentAnalyticField.cc.

240 : ComponentBase() {
241 m_className = "ComponentAnalyticField";
242 CellInit();
243}
Garfield::ComponentBase::ComponentBase
ComponentBase()
Constructor.
Definition: ComponentBase.cc:9
Garfield::ComponentBase::m_className
std::string m_className
Class name.
Definition: ComponentBase.hh:285

◆ ~ComponentAnalyticField()

Garfield::ComponentAnalyticField::~ComponentAnalyticField ( )
inline

Destructor.

Definition at line 20 of file ComponentAnalyticField.hh.

20{}

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 1581 of file ComponentAnalyticField.cc.

1582 {
1583 // Convert from fC to internal units (division by 4 pi epsilon0).
1584 Charge3d charge;
1585 charge.x = x;
1586 charge.y = y;
1587 charge.z = z;
1588 charge.e = q / FourPiEpsilon0;
1589 m_ch3d.push_back(std::move(charge));
1590}

◆ 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 1373 of file ComponentAnalyticField.cc.

1376 {
1377 if (!m_polar || (!m_ynplan[2] && !m_ynplan[3])) {
1378 std::cerr << m_className << "::AddPixelOnPlanePhi:\n"
1379 << " There are no planes at constant phi.\n";
1380 return;
1381 }
1382
1383 if (fabs(rmax - rmin) < Small || fabs(zmax - zmin) < Small) {
1384 std::cerr << m_className << "::AddPixelOnPlanePhi:\n"
1385 << " Pixel width must be greater than zero.\n";
1386 return;
1387 }
1388 if (rmin < Small || rmax < Small) {
1389 std::cerr << m_className << "::AddPixelOnPlanePhi:\n"
1390 << " Radius must be greater than zero.\n";
1391 return;
1392 }
1393 Pixel newPixel;
1394 newPixel.type = label;
1395 newPixel.ind = -1;
1396 const double smin = log(rmin);
1397 const double smax = log(rmax);
1398 newPixel.smin = std::min(smin, smax);
1399 newPixel.smax = std::max(smin, smax);
1400 newPixel.zmin = std::min(zmin, zmax);
1401 newPixel.zmax = std::max(zmin, zmax);
1402 newPixel.gap = gap > Small ? DegreeToRad * gap : -1.;
1403
1404 int iplane = 2;
1405 if (m_ynplan[3]) {
1406 const double d0 = fabs(m_coplan[2] - phi * DegreeToRad);
1407 const double d1 = fabs(m_coplan[3] - phi * DegreeToRad);
1408 if (d1 < d0) iplane = 3;
1409 }
1410
1411 m_planes[iplane].pixels.push_back(std::move(newPixel));
1412}
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 1335 of file ComponentAnalyticField.cc.

1338 {
1339 if (!m_polar || (!m_ynplan[0] && !m_ynplan[1])) {
1340 std::cerr << m_className << "::AddPixelOnPlaneR:\n"
1341 << " There are no planes at constant r.\n";
1342 return;
1343 }
1344
1345 if (fabs(phimax - phimin) < Small || fabs(zmax - zmin) < Small) {
1346 std::cerr << m_className << "::AddPixelOnPlaneR:\n"
1347 << " Pixel width must be greater than zero.\n";
1348 return;
1349 }
1350
1351 Pixel newPixel;
1352 newPixel.type = label;
1353 newPixel.ind = -1;
1354 const double smin = phimin * DegreeToRad;
1355 const double smax = phimax * DegreeToRad;
1356 newPixel.smin = std::min(smin, smax);
1357 newPixel.smax = std::max(smin, smax);
1358 newPixel.zmin = std::min(zmin, zmax);
1359 newPixel.zmax = std::max(zmin, zmax);
1360 newPixel.gap = gap > Small ? gap : -1.;
1361
1362 int iplane = 0;
1363 if (m_ynplan[1]) {
1364 const double rho = r > 0. ? log(r) : -25.;
1365 const double d0 = fabs(m_coplan[0] - rho);
1366 const double d1 = fabs(m_coplan[1] - rho);
1367 if (d1 < d0) iplane = 1;
1368 }
1369
1370 m_planes[iplane].pixels.push_back(std::move(newPixel));
1371}

◆ 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. 
)

Add a pixel on an existing plane at constant x.

Definition at line 1267 of file ComponentAnalyticField.cc.

1269 {
1270 if (m_polar || (!m_ynplan[0] && !m_ynplan[1])) {
1271 std::cerr << m_className << "::AddPixelOnPlaneX:\n"
1272 << " There are no planes at constant x.\n";
1273 return;
1274 }
1275
1276 if (fabs(ymax - ymin) < Small || fabs(zmax - zmin) < Small) {
1277 std::cerr << m_className << "::AddPixelOnPlaneX:\n"
1278 << " Pixel width must be greater than zero.\n";
1279 return;
1280 }
1281
1282 Pixel newPixel;
1283 newPixel.type = label;
1284 newPixel.ind = -1;
1285 newPixel.smin = std::min(ymin, ymax);
1286 newPixel.smax = std::max(ymin, ymax);
1287 newPixel.zmin = std::min(zmin, zmax);
1288 newPixel.zmax = std::max(zmin, zmax);
1289 newPixel.gap = gap > Small ? gap : -1.;
1290
1291 int iplane = 0;
1292 if (m_ynplan[1]) {
1293 const double d0 = fabs(m_coplan[0] - x);
1294 const double d1 = fabs(m_coplan[1] - x);
1295 if (d1 < d0) iplane = 1;
1296 }
1297
1298 m_planes[iplane].pixels.push_back(std::move(newPixel));
1299}

◆ 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. 
)

Add a pixel on an existing plane at constant y.

Definition at line 1301 of file ComponentAnalyticField.cc.

1303 {
1304 if (m_polar || (!m_ynplan[2] && !m_ynplan[3])) {
1305 std::cerr << m_className << "::AddPixelOnPlaneY:\n"
1306 << " There are no planes at constant y.\n";
1307 return;
1308 }
1309
1310 if (fabs(xmax - xmin) < Small || fabs(zmax - zmin) < Small) {
1311 std::cerr << m_className << "::AddPixelOnPlaneY:\n"
1312 << " Pixel width must be greater than zero.\n";
1313 return;
1314 }
1315
1316 Pixel newPixel;
1317 newPixel.type = label;
1318 newPixel.ind = -1;
1319 newPixel.smin = std::min(xmin, xmax);
1320 newPixel.smax = std::max(xmin, xmax);
1321 newPixel.zmin = std::min(zmin, zmax);
1322 newPixel.zmax = std::max(zmin, zmax);
1323 newPixel.gap = gap > Small ? gap : -1.;
1324
1325 int iplane = 2;
1326 if (m_ynplan[3]) {
1327 const double d0 = fabs(m_coplan[2] - y);
1328 const double d1 = fabs(m_coplan[3] - y);
1329 if (d1 < d0) iplane = 3;
1330 }
1331
1332 m_planes[iplane].pixels.push_back(std::move(newPixel));
1333}

◆ AddPlanePhi()

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

Add a plane at constant phi.

Definition at line 1029 of file ComponentAnalyticField.cc.

1030 {
1031 if (!m_polar) {
1032 std::cerr << m_className << "::AddPlanePhi:\n"
1033 << " Not compatible with Cartesian coordinates; ignored.\n";
1034 return;
1035 }
1036 if (m_ynplan[2] && m_ynplan[3]) {
1037 std::cerr << m_className << "::AddPlanePhi:\n"
1038 << " Cannot have more than two planes at constant phi.\n";
1039 return;
1040 }
1041
1042 if (m_ynplan[2]) {
1043 m_ynplan[3] = true;
1044 m_coplan[3] = phi * DegreeToRad;
1045 m_vtplan[3] = v;
1046 m_planes[3].type = label;
1047 m_planes[3].ind = -1;
1048 } else {
1049 m_ynplan[2] = true;
1050 m_coplan[2] = phi * DegreeToRad;
1051 m_vtplan[2] = v;
1052 m_planes[2].type = label;
1053 m_planes[2].ind = -1;
1054 // Switch off default periodicity.
1055 if (m_pery && std::abs(m_sy - TwoPi) < 1.e-4) {
1056 m_pery = false;
1057 }
1058 }
1059
1060 // Force recalculation of the capacitance and signal matrices.
1061 m_cellset = false;
1062 m_sigset = false;
1063}

◆ AddPlaneR()

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

Add a plane at constant radius.

Definition at line 991 of file ComponentAnalyticField.cc.

992 {
993 if (!m_polar) {
994 std::cerr << m_className << "::AddPlaneR:\n"
995 << " Not compatible with Cartesian coordinates; ignored.\n";
996 return;
997 }
998 if (r <= 0.) {
999 std::cerr << m_className << "::AddPlaneR:\n"
1000 << " Radius must be larger than zero; plane ignored.\n";
1001 return;
1002 }
1003
1004 if (m_ynplan[0] && m_ynplan[1]) {
1005 std::cerr << m_className << "::AddPlaneR:\n"
1006 << " Cannot have more than two circular planes.\n";
1007 return;
1008 }
1009
1010 if (m_ynplan[0]) {
1011 m_ynplan[1] = true;
1012 m_coplan[1] = log(r);
1013 m_vtplan[1] = v;
1014 m_planes[1].type = label;
1015 m_planes[1].ind = -1;
1016 } else {
1017 m_ynplan[0] = true;
1018 m_coplan[0] = log(r);
1019 m_vtplan[0] = v;
1020 m_planes[0].type = label;
1021 m_planes[0].ind = -1;
1022 }
1023
1024 // Force recalculation of the capacitance and signal matrices.
1025 m_cellset = false;
1026 m_sigset = false;
1027}

◆ AddPlaneX()

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

Add a plane at constant x.

Definition at line 927 of file ComponentAnalyticField.cc.

928 {
929 if (m_polar) {
930 std::cerr << m_className << "::AddPlaneX:\n"
931 << " Not compatible with polar coordinates; ignored.\n";
932 return;
933 }
934 if (m_ynplan[0] && m_ynplan[1]) {
935 std::cerr << m_className << "::AddPlaneX:\n"
936 << " Cannot have more than two planes at constant x.\n";
937 return;
938 }
939
940 if (m_ynplan[0]) {
941 m_ynplan[1] = true;
942 m_coplan[1] = x;
943 m_vtplan[1] = v;
944 m_planes[1].type = label;
945 m_planes[1].ind = -1;
946 } else {
947 m_ynplan[0] = true;
948 m_coplan[0] = x;
949 m_vtplan[0] = v;
950 m_planes[0].type = label;
951 m_planes[0].ind = -1;
952 }
953
954 // Force recalculation of the capacitance and signal matrices.
955 m_cellset = false;
956 m_sigset = false;
957}

◆ AddPlaneY()

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

Add a plane at constant y.

Definition at line 959 of file ComponentAnalyticField.cc.

960 {
961 if (m_polar) {
962 std::cerr << m_className << "::AddPlaneY:\n"
963 << " Not compatible with polar coordinates; ignored.\n";
964 return;
965 }
966 if (m_ynplan[2] && m_ynplan[3]) {
967 std::cerr << m_className << "::AddPlaneY:\n"
968 << " Cannot have more than two planes at constant y.\n";
969 return;
970 }
971
972 if (m_ynplan[2]) {
973 m_ynplan[3] = true;
974 m_coplan[3] = y;
975 m_vtplan[3] = v;
976 m_planes[3].type = label;
977 m_planes[3].ind = -1;
978 } else {
979 m_ynplan[2] = true;
980 m_coplan[2] = y;
981 m_vtplan[2] = v;
982 m_planes[2].type = label;
983 m_planes[2].ind = -1;
984 }
985
986 // Force recalculation of the capacitance and signal matrices.
987 m_cellset = false;
988 m_sigset = false;
989}

◆ AddReadout()

void Garfield::ComponentAnalyticField::AddReadout ( const std::string &  label)

Setup the weighting field for a given group of wires or planes.

Definition at line 3512 of file ComponentAnalyticField.cc.

3512 {
3513 // Check if this readout group already exists.
3514 if (std::find(m_readout.begin(), m_readout.end(), label) != m_readout.end()) {
3515 std::cout << m_className << "::AddReadout:\n";
3516 std::cout << " Readout group " << label << " already exists.\n";
3517 return;
3518 }
3519 m_readout.push_back(label);
3520
3521 unsigned int nWiresFound = 0;
3522 for (const auto& wire : m_w) {
3523 if (wire.type == label) ++nWiresFound;
3524 }
3525
3526 unsigned int nPlanesFound = 0;
3527 unsigned int nStripsFound = 0;
3528 unsigned int nPixelsFound = 0;
3529 for (int i = 0; i < 5; ++i) {
3530 if (m_planes[i].type == label) ++nPlanesFound;
3531 for (const auto& strip : m_planes[i].strips1) {
3532 if (strip.type == label) ++nStripsFound;
3533 }
3534 for (const auto& strip : m_planes[i].strips2) {
3535 if (strip.type == label) ++nStripsFound;
3536 }
3537 for (const auto& pixel : m_planes[i].pixels) {
3538 if (pixel.type == label) ++nPixelsFound;
3539 }
3540 }
3541
3542 if (nWiresFound == 0 && nPlanesFound == 0 && nStripsFound == 0 &&
3543 nPixelsFound == 0) {
3544 std::cerr << m_className << "::AddReadout:\n";
3545 std::cerr << " At present there are no wires, planes or strips\n";
3546 std::cerr << " associated to readout group " << label << ".\n";
3547 } else {
3548 std::cout << m_className << "::AddReadout:\n";
3549 std::cout << " Readout group " << label << " comprises:\n";
3550 if (nWiresFound > 1) {
3551 std::cout << " " << nWiresFound << " wires\n";
3552 } else if (nWiresFound == 1) {
3553 std::cout << " 1 wire\n";
3554 }
3555 if (nPlanesFound > 1) {
3556 std::cout << " " << nPlanesFound << " planes\n";
3557 } else if (nPlanesFound == 1) {
3558 std::cout << " 1 plane\n";
3559 }
3560 if (nStripsFound > 1) {
3561 std::cout << " " << nStripsFound << " strips\n";
3562 } else if (nStripsFound == 1) {
3563 std::cout << " 1 strip\n";
3564 }
3565 if (nPixelsFound > 1) {
3566 std::cout << " " << nPixelsFound << " pixels\n";
3567 } else if (nPixelsFound == 1) {
3568 std::cout << " 1 pixel\n";
3569 }
3570 }
3571
3572 m_sigset = false;
3573}

◆ 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 1208 of file ComponentAnalyticField.cc.

1213 {
1214 if (!m_polar || (!m_ynplan[2] && !m_ynplan[3])) {
1215 std::cerr << m_className << "::AddStripOnPlanePhi:\n"
1216 << " There are no planes at constant phi.\n";
1217 return;
1218 }
1219
1220 if (dir != 'r' && dir != 'R' && dir != 'z' && dir != 'Z') {
1221 std::cerr << m_className << "::AddStripOnPlanePhi:\n"
1222 << " Invalid direction (" << dir << ").\n"
1223 << " Only strips in r or z direction are possible.\n";
1224 return;
1225 }
1226
1227 if (fabs(smax - smin) < Small) {
1228 std::cerr << m_className << "::AddStripOnPlanePhi:\n"
1229 << " Strip width must be greater than zero.\n";
1230 return;
1231 }
1232
1233 Strip newStrip;
1234 newStrip.type = label;
1235 newStrip.ind = -1;
1236 if (dir== 'z' || dir == 'Z') {
1237 if (smin < Small || smax < Small) {
1238 std::cerr << m_className << "::AddStripOnPlanePhi:\n"
1239 << " Radius must be greater than zero.\n";
1240 return;
1241 }
1242 const double rhomin = log(smin);
1243 const double rhomax = log(smax);
1244 newStrip.smin = std::min(rhomin, rhomax);
1245 newStrip.smax = std::max(rhomin, rhomax);
1246 } else {
1247 newStrip.smin = std::min(smin, smax);
1248 newStrip.smax = std::max(smin, smax);
1249 }
1250 newStrip.gap = gap > Small ? DegreeToRad * gap : -1.;
1251
1252 int iplane = 2;
1253 if (m_ynplan[3]) {
1254 const double d2 = fabs(m_coplan[2] - phi * DegreeToRad);
1255 const double d3 = fabs(m_coplan[3] - phi * DegreeToRad);
1256 if (d3 < d2) iplane = 3;
1257 }
1258
1259 if (dir == 'r' || dir == 'R') {
1260 m_planes[iplane].strips1.push_back(std::move(newStrip));
1261 } else {
1262 m_planes[iplane].strips2.push_back(std::move(newStrip));
1263 }
1264}

◆ 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 1155 of file ComponentAnalyticField.cc.

1159 {
1160 if (!m_polar || (!m_ynplan[0] && !m_ynplan[1])) {
1161 std::cerr << m_className << "::AddStripOnPlaneR:\n"
1162 << " There are no planes at constant r.\n";
1163 return;
1164 }
1165
1166 if (dir != 'p' && dir != 'P' && dir != 'z' && dir != 'Z') {
1167 std::cerr << m_className << "::AddStripOnPlaneR:\n"
1168 << " Invalid direction (" << dir << ").\n"
1169 << " Only strips in p(hi) or z direction are possible.\n";
1170 return;
1171 }
1172
1173 if (fabs(smax - smin) < Small) {
1174 std::cerr << m_className << "::AddStripOnPlaneR:\n"
1175 << " Strip width must be greater than zero.\n";
1176 return;
1177 }
1178
1179 Strip newStrip;
1180 newStrip.type = label;
1181 newStrip.ind = -1;
1182 if (dir == 'z' || dir == 'Z') {
1183 const double phimin = smin * DegreeToRad;
1184 const double phimax = smax * DegreeToRad;
1185 newStrip.smin = std::min(phimin, phimax);
1186 newStrip.smax = std::max(phimin, phimax);
1187 } else {
1188 newStrip.smin = std::min(smin, smax);
1189 newStrip.smax = std::max(smin, smax);
1190 }
1191 newStrip.gap = gap > Small ? gap : -1.;
1192
1193 int iplane = 0;
1194 if (m_ynplan[1]) {
1195 const double rho = r > 0. ? log(r) : -25.;
1196 const double d0 = fabs(m_coplan[0] - rho);
1197 const double d1 = fabs(m_coplan[1] - rho);
1198 if (d1 < d0) iplane = 1;
1199 }
1200
1201 if (dir == 'p' || dir == 'P') {
1202 m_planes[iplane].strips1.push_back(std::move(newStrip));
1203 } else {
1204 m_planes[iplane].strips2.push_back(std::move(newStrip));
1205 }
1206}

◆ 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.

Definition at line 1065 of file ComponentAnalyticField.cc.

1069 {
1070 if (m_polar || (!m_ynplan[0] && !m_ynplan[1])) {
1071 std::cerr << m_className << "::AddStripOnPlaneX:\n"
1072 << " There are no planes at constant x.\n";
1073 return;
1074 }
1075
1076 if (dir != 'y' && dir != 'Y' && dir != 'z' && dir != 'Z') {
1077 std::cerr << m_className << "::AddStripOnPlaneX:\n"
1078 << " Invalid direction (" << dir << ").\n"
1079 << " Only strips in y or z direction are possible.\n";
1080 return;
1081 }
1082
1083 if (fabs(smax - smin) < Small) {
1084 std::cerr << m_className << "::AddStripOnPlaneX:\n"
1085 << " Strip width must be greater than zero.\n";
1086 return;
1087 }
1088
1089 Strip newStrip;
1090 newStrip.type = label;
1091 newStrip.ind = -1;
1092 newStrip.smin = std::min(smin, smax);
1093 newStrip.smax = std::max(smin, smax);
1094 newStrip.gap = gap > Small ? gap : -1.;
1095
1096 int iplane = 0;
1097 if (m_ynplan[1]) {
1098 const double d0 = fabs(m_coplan[0] - x);
1099 const double d1 = fabs(m_coplan[1] - x);
1100 if (d1 < d0) iplane = 1;
1101 }
1102
1103 if (dir == 'y' || dir == 'Y') {
1104 m_planes[iplane].strips1.push_back(std::move(newStrip));
1105 } else {
1106 m_planes[iplane].strips2.push_back(std::move(newStrip));
1107 }
1108}

◆ 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 1110 of file ComponentAnalyticField.cc.

1114 {
1115 if (m_polar || (!m_ynplan[2] && !m_ynplan[3])) {
1116 std::cerr << m_className << "::AddStripOnPlaneY:\n"
1117 << " There are no planes at constant y.\n";
1118 return;
1119 }
1120
1121 if (dir != 'x' && dir != 'X' && dir != 'z' && dir != 'Z') {
1122 std::cerr << m_className << "::AddStripOnPlaneY:\n"
1123 << " Invalid direction (" << dir << ").\n"
1124 << " Only strips in x or z direction are possible.\n";
1125 return;
1126 }
1127
1128 if (fabs(smax - smin) < Small) {
1129 std::cerr << m_className << "::AddStripOnPlaneY:\n"
1130 << " Strip width must be greater than zero.\n";
1131 return;
1132 }
1133
1134 Strip newStrip;
1135 newStrip.type = label;
1136 newStrip.ind = -1;
1137 newStrip.smin = std::min(smin, smax);
1138 newStrip.smax = std::max(smin, smax);
1139 newStrip.gap = gap > Small ? gap : -1.;
1140
1141 int iplane = 2;
1142 if (m_ynplan[3]) {
1143 const double d2 = fabs(m_coplan[2] - y);
1144 const double d3 = fabs(m_coplan[3] - y);
1145 if (d3 < d2) iplane = 3;
1146 }
1147
1148 if (dir == 'x' || dir == 'X') {
1149 m_planes[iplane].strips1.push_back(std::move(newStrip));
1150 } else {
1151 m_planes[iplane].strips2.push_back(std::move(newStrip));
1152 }
1153}

◆ AddTube()

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

Add a tube.

Definition at line 887 of file ComponentAnalyticField.cc.

889 {
890 // Check if the provided parameters make sense.
891 if (radius <= 0.0) {
892 std::cerr << m_className << "::AddTube: Unphysical tube dimension.\n";
893 return;
894 }
895
896 if (nEdges < 3 && nEdges != 0) {
897 std::cerr << m_className << "::AddTube: Unphysical number of tube edges ("
898 << nEdges << ")\n";
899 return;
900 }
901
902 // If there is already a tube defined, print a warning message.
903 if (m_tube) {
904 std::cout << m_className << "::AddTube:\n"
905 << " Warning: Existing tube settings will be overwritten.\n";
906 }
907
908 // Set the coordinate system.
909 m_tube = true;
910 m_polar = false;
911
912 // Set the tube parameters.
913 m_cotube = radius;
914 m_cotube2 = radius * radius;
915 m_vttube = voltage;
916
917 m_ntube = nEdges;
918
919 m_planes[4].type = label;
920 m_planes[4].ind = -1;
921
922 // Force recalculation of the capacitance and signal matrices.
923 m_cellset = false;
924 m_sigset = false;
925}

Referenced by GarfieldPhysics::CreateGeometry().

◆ 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 820 of file ComponentAnalyticField.cc.

825 {
826 // Check if the provided parameters make sense.
827 if (diameter <= 0.) {
828 std::cerr << m_className << "::AddWire: Unphysical wire diameter.\n";
829 return;
830 }
831
832 if (tension <= 0.) {
833 std::cerr << m_className << "::AddWire: Unphysical wire tension.\n";
834 return;
835 }
836
837 if (rho <= 0.) {
838 std::cerr << m_className << "::AddWire: Unphysical wire density.\n";
839 return;
840 }
841
842 if (length <= 0.) {
843 std::cerr << m_className << "::AddWire: Unphysical wire length.\n";
844 return;
845 }
846
847 if (ntrap <= 0) {
848 std::cerr << m_className << "::AddWire: Nbr. of trap radii must be > 0.\n";
849 return;
850 }
851
852 if (m_polar && x <= diameter) {
853 std::cerr << m_className << "::AddWire: Wire is too close to the origin.\n";
854 return;
855 }
856
857 // Create a new wire.
858 Wire wire;
859 if (m_polar) {
860 double r = 0., phi = 0.;
861 Polar2Internal(x, y, r, phi);
862 wire.x = r;
863 wire.y = phi;
864 wire.r = 0.5 * diameter / x;
865 } else {
866 wire.x = x;
867 wire.y = y;
868 wire.r = 0.5 * diameter;
869 }
870 wire.v = voltage;
871 wire.u = length;
872 wire.type = label;
873 wire.e = 0.;
874 wire.ind = -1;
875 wire.nTrap = ntrap;
876 wire.tension = tension;
877 wire.density = rho;
878 // Add the wire to the list
879 m_w.push_back(std::move(wire));
880 ++m_nWires;
881
882 // Force recalculation of the capacitance and signal matrices.
883 m_cellset = false;
884 m_sigset = false;
885}

Referenced by GarfieldPhysics::CreateGeometry().

◆ ClearCharges()

void Garfield::ComponentAnalyticField::ClearCharges ( )

Remove all point charges.

Definition at line 1592 of file ComponentAnalyticField.cc.

1592 {
1593 m_ch3d.clear();
1594 m_nTermBessel = 10;
1595 m_nTermPoly = 100;
1596}

◆ ElectricField() [1/2]

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::ComponentBase.

Definition at line 40 of file ComponentAnalyticField.hh.

42 {
43 m = nullptr;
44 // Calculate the field.
45 status = Field(x, y, z, ex, ey, ez, v, true);
46 // If the field is ok, get the medium.
47 if (status == 0) {
48 m = m_geometry ? m_geometry->GetMedium(x, y, z) : m_medium;
49 if (!m) {
50 status = -6;
51 } else if (!m->IsDriftable()) {
52 status = -5;
53 }
54 }
55 }
Garfield::ComponentBase::m_geometry
GeometryBase * m_geometry
Pointer to the geometry.
Definition: ComponentBase.hh:288
Garfield::GeometryBase::GetMedium
virtual Medium * GetMedium(const double x, const double y, const double z) const =0
Retrieve the medium at a given point.

◆ ElectricField() [2/2]

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::ComponentBase.

Definition at line 23 of file ComponentAnalyticField.hh.

24 {
25 m = nullptr;
26 // Calculate the field.
27 double v = 0.;
28 status = Field(x, y, z, ex, ey, ez, v, false);
29 // If the field is ok, get the medium.
30 if (status == 0) {
31 m = m_geometry ? m_geometry->GetMedium(x, y, z) : m_medium;
32 if (!m) {
33 status = -6;
34 } else if (!m->IsDriftable()) {
35 status = -5;
36 }
37 }
38 }

◆ 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 1750 of file ComponentAnalyticField.cc.

1751 {
1752 //-----------------------------------------------------------------------
1753 // FFIELD - Subroutine calculating the electric field at a given wire
1754 // position, as if the wire itself were not there but with
1755 // the presence of its mirror images.
1756 // VARIABLES : IW : wire number
1757 // EX, EY : x- and y-component of the electric field.
1758 // (Last changed on 27/ 1/96.)
1759 //-----------------------------------------------------------------------
1760 ex = ey = 0.;
1761 // Check the wire number.
1762 if (iw >= m_nWires) {
1763 std::cerr << m_className << "::ElectricFieldAtWire: Index out of range.\n";
1764 return false;
1765 }
1766 // Set the flags appropriately.
1767 std::vector<bool> cnalso(m_nWires, true);
1768 cnalso[iw] = false;
1769
1770 const double xpos = m_w[iw].x;
1771 const double ypos = m_w[iw].y;
1772 // Call the appropriate function.
1773 switch (m_cellType) {
1774 case A00:
1775 FieldAtWireA00(xpos, ypos, ex, ey, cnalso);
1776 break;
1777 case B1X:
1778 FieldAtWireB1X(xpos, ypos, ex, ey, cnalso);
1779 break;
1780 case B1Y:
1781 FieldAtWireB1Y(xpos, ypos, ex, ey, cnalso);
1782 break;
1783 case B2X:
1784 FieldAtWireB2X(xpos, ypos, ex, ey, cnalso);
1785 break;
1786 case B2Y:
1787 FieldAtWireB2Y(xpos, ypos, ex, ey, cnalso);
1788 break;
1789 case C10:
1790 FieldAtWireC10(xpos, ypos, ex, ey, cnalso);
1791 break;
1792 case C2X:
1793 FieldAtWireC2X(xpos, ypos, ex, ey, cnalso);
1794 break;
1795 case C2Y:
1796 FieldAtWireC2Y(xpos, ypos, ex, ey, cnalso);
1797 break;
1798 case C30:
1799 FieldAtWireC30(xpos, ypos, ex, ey, cnalso);
1800 break;
1801 case D10:
1802 FieldAtWireD10(xpos, ypos, ex, ey, cnalso);
1803 break;
1804 case D20:
1805 FieldAtWireD20(xpos, ypos, ex, ey, cnalso);
1806 break;
1807 case D30:
1808 FieldAtWireD30(xpos, ypos, ex, ey, cnalso);
1809 break;
1810 default:
1811 std::cerr << m_className << "::ElectricFieldAtWire:\n"
1812 << " Unknown cell type (id " << m_cellType << ")\n";
1813 return false;
1814 break;
1815 }
1816 // Correct for the equipotential planes.
1817 ex -= m_corvta;
1818 ey -= m_corvtb;
1819 if (m_polar) {
1820 const double r = exp(xpos);
1821 const double er = ex / r;
1822 const double ep = ey / r;
1823 const double ct = cos(ypos);
1824 const double st = sin(ypos);
1825 ex = +ct * er - st * ep;
1826 ey = +st * er + ct * ep;
1827 }
1828 return true;
1829}
Heed::cos
DoubleAc cos(const DoubleAc &f)
Definition: DoubleAc.cpp:432
Heed::exp
DoubleAc exp(const DoubleAc &f)
Definition: DoubleAc.cpp:377
Heed::sin
DoubleAc sin(const DoubleAc &f)
Definition: DoubleAc.cpp:384

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 223 of file ComponentAnalyticField.hh.

223{ 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 1414 of file ComponentAnalyticField.cc.

1414 {
1415
1416 m_cellset = false;
1417 m_sigset = false;
1418 m_dipole = on;
1419}

◆ 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 275 of file ComponentAnalyticField.hh.

275{ m_extrapolateForces = 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 1895 of file ComponentAnalyticField.cc.

1898 {
1899 if (!m_cellset && !Prepare()) {
1900 std::cerr << m_className << "::ForcesOnWire: Cell not set up.\n";
1901 return false;
1902 }
1903
1904 if (iw >= m_nWires) {
1905 std::cerr << m_className << "::ForcesOnWire: Wire index out of range.\n";
1906 return false;
1907 }
1908 if (m_polar) {
1909 std::cerr << m_className << "::ForcesOnWire: Cannot handle polar cells.\n";
1910 return false;
1911 }
1912 const auto& wire = m_w[iw];
1913 // Compute a 'safe box' around the wire.
1914 double bxmin = m_perx ? wire.x - 0.5 * m_sx : 2 * m_xmin - m_xmax;
1915 double bxmax = m_perx ? wire.x + 0.5 * m_sx : 2 * m_xmax - m_xmin;
1916 double bymin = m_pery ? wire.y - 0.5 * m_sy : 2 * m_ymin - m_ymax;
1917 double bymax = m_pery ? wire.y + 0.5 * m_sy : 2 * m_ymax - m_ymin;
1918
1919 // If the initial area is almost zero in 1 direction, make it square.
1920 if (std::abs(bxmax - bxmin) < 0.1 * std::abs(bymax - bymin)) {
1921 bxmin = wire.x - 0.5 * std::abs(bymax - bymin);
1922 bxmax = wire.x + 0.5 * std::abs(bymax - bymin);
1923 } else if (std::abs(bymax - bymin) < 0.1 * std::abs(bxmax - bxmin)) {
1924 bymin = wire.y - 0.5 * std::abs(bxmax - bxmin);
1925 bymax = wire.y + 0.5 * std::abs(bxmax - bxmin);
1926 }
1927 const double dw = 2 * wire.r;
1928 // Scan the other wires.
1929 for (unsigned int j = 0; j < m_nWires; ++j) {
1930 if (j == iw) continue;
1931 const double xj = m_w[j].x;
1932 const double yj = m_w[j].y;
1933 const double dj = 2 * m_w[j].r;
1934 double xnear = m_perx ? xj - m_sx * int(round((xj - wire.x) / m_sx)) : xj;
1935 double ynear = m_pery ? yj - m_sy * int(round((yj - wire.y) / m_sy)) : yj;
1936 if (std::abs(xnear - wire.x) > std::abs(ynear - wire.y)) {
1937 if (xnear < wire.x) {
1938 bxmin = std::max(bxmin, xnear + dj + dw);
1939 if (m_perx) bxmax = std::min(bxmax, xnear + m_sx - dj - dw);
1940 } else {
1941 bxmax = std::min(bxmax, xnear - dj - dw);
1942 if (m_perx) bxmin = std::max(bxmin, xnear - m_sx + dj + dw);
1943 }
1944 } else {
1945 if (ynear < wire.y) {
1946 bymin = std::max({bymin, ynear - dj - dw, ynear + dj + dw});
1947 if (m_pery) bymax = std::min(bymax, ynear + m_sy - dj - dw);
1948 } else {
1949 bymax = std::min({bymax, ynear - dj - dw, ynear + dj + dw});
1950 if (m_pery) bymin = std::max(bymin, ynear - m_sy + dj + dw);
1951 }
1952 }
1953 }
1954 // Scan the planes.
1955 if (m_ynplan[0]) bxmin = std::max(bxmin, m_coplan[0] + dw);
1956 if (m_ynplan[1]) bxmax = std::min(bxmax, m_coplan[1] - dw);
1957 if (m_ynplan[2]) bymin = std::max(bymin, m_coplan[2] + dw);
1958 if (m_ynplan[3]) bymax = std::min(bymax, m_coplan[3] - dw);
1959
1960 // If there is a tube, check all corners.
1961 if (m_tube) {
1962 const double d2 = m_cotube2 - dw * dw;
1963 if (d2 < Small) {
1964 std::cerr << m_className << "::ForcesOnWire:\n Diameter of wire " << iw
1965 << " is too large compared to the tube.\n";
1966 return false;
1967 }
1968
1969 double corr = sqrt((bxmin * bxmin + bymin * bymin) / d2);
1970 if (corr > 1.) {
1971 bxmin /= corr;
1972 bymin /= corr;
1973 }
1974 corr = sqrt((bxmin * bxmin + bymax * bymax) / d2);
1975 if (corr > 1.) {
1976 bxmin /= corr;
1977 bymax /= corr;
1978 }
1979 corr = sqrt((bxmax * bxmax + bymin * bymin) / d2);
1980 if (corr > 1.) {
1981 bxmax /= corr;
1982 bymin /= corr;
1983 }
1984 corr = sqrt((bxmax * bxmax + bymax * bymax) / d2);
1985 if (corr > 1) {
1986 bxmax /= corr;
1987 bymax /= corr;
1988 }
1989 }
1990 // Make sure we found a reasonable 'safe area'.
1991 if ((bxmin - wire.x) * (wire.x - bxmax) <= 0 ||
1992 (bymin - wire.y) * (wire.y - bymax) <= 0) {
1993 std::cerr << m_className << "::ForcesOnWire:\n Unable to find an area "
1994 << "free of elements around wire " << iw << ".\n";
1995 return false;
1996 }
1997 // Now set a reasonable scanning range.
1998 double sxmin = bxmin;
1999 double sxmax = bxmax;
2000 double symin = bymin;
2001 double symax = bymax;
2002 if (m_scanRange == ScanningRange::User) {
2003 // User-specified range.
2004 sxmin = wire.x + m_xScanMin;
2005 symin = wire.y + m_yScanMin;
2006 sxmax = wire.x + m_xScanMax;
2007 symax = wire.y + m_yScanMax;
2008 } else if (m_scanRange == ScanningRange::FirstOrder) {
2009 // Get the field and force at the nominal position.
2010 double ex = 0., ey = 0.;
2011 ElectricFieldAtWire(iw, ex, ey);
2012 double fx = -ex * wire.e * Internal2Newton;
2013 double fy = -ey * wire.e * Internal2Newton;
2014 if (m_useGravitationalForce) {
2015 // Mass per unit length [kg / cm].
2016 const double m = 1.e-3 * wire.density * Pi * wire.r * wire.r;
2017 fx -= m_down[0] * GravitationalAcceleration * m;
2018 fy -= m_down[1] * GravitationalAcceleration * m;
2019 }
2020 const double u2 = wire.u * wire.u;
2021 const double shiftx =
2022 -125 * fx * u2 / (GravitationalAcceleration * wire.tension);
2023 const double shifty =
2024 -125 * fy * u2 / (GravitationalAcceleration * wire.tension);
2025 // If 0th order estimate of shift is not small.
2026 const double tol = 0.1 * wire.r;
2027 if (std::abs(shiftx) > tol || std::abs(shifty) > tol) {
2028 sxmin = std::max(bxmin, std::min(wire.x + m_scaleRange * shiftx,
2029 wire.x - shiftx / m_scaleRange));
2030 symin = std::max(bymin, std::min(wire.y + m_scaleRange * shifty,
2031 wire.y - shifty / m_scaleRange));
2032 sxmax = std::min(bxmax, std::max(wire.x + m_scaleRange * shiftx,
2033 wire.x - shiftx / m_scaleRange));
2034 symax = std::min(bymax, std::max(wire.y + m_scaleRange * shifty,
2035 wire.y - shifty / m_scaleRange));
2036 // If one is very small, make the area square within bounds.
2037 if (std::abs(sxmax - sxmin) < 0.1 * std::abs(symax - symin)) {
2038 sxmin = std::max(bxmin, wire.x - 0.5 * std::abs(symax - symin));
2039 sxmax = std::min(bxmax, wire.x + 0.5 * std::abs(symax - symin));
2040 } else if (std::abs(symax - symin) < 0.1 * std::abs(sxmax - sxmin)) {
2041 symin = std::max(bymin, wire.y - 0.5 * std::abs(sxmax - sxmin));
2042 symax = std::min(bymax, wire.y + 0.5 * std::abs(sxmax - sxmin));
2043 }
2044 }
2045 }
2046 if (m_debug) {
2047 std::cout << m_className << "::ForcesOnWire:\n";
2048 std::printf(" Free area %12.5e < x < %12.5e\n", bxmin, bxmax);
2049 std::printf(" %12.5e < y < %12.5e\n", bymin, bymax);
2050 std::printf(" Scan area %12.5e < x < %12.5e\n", sxmin, sxmax);
2051 std::printf(" %12.5e < y < %12.5e\n", symin, symax);
2052 }
2053
2054 xMap.resize(m_nScanX);
2055 const double stepx = (sxmax - sxmin) / (m_nScanX - 1);
2056 for (unsigned int i = 0; i < m_nScanX; ++i) {
2057 xMap[i] = sxmin + i * stepx;
2058 }
2059 yMap.resize(m_nScanY);
2060 const double stepy = (symax - symin) / (m_nScanY - 1);
2061 for (unsigned int i = 0; i < m_nScanY; ++i) {
2062 yMap[i] = symin + i * stepy;
2063 }
2064 // Save the original coordinates of the wire.
2065 const double x0 = wire.x;
2066 const double y0 = wire.y;
2067 // Prepare interpolation tables.
2068 fxMap.assign(m_nScanX, std::vector<double>(m_nScanY, 0.));
2069 fyMap.assign(m_nScanX, std::vector<double>(m_nScanY, 0.));
2070 bool ok = true;
2071 for (unsigned int i = 0; i < m_nScanX; ++i) {
2072 for (unsigned int j = 0; j < m_nScanY; ++j) {
2073 // Get the wire position for this shift.
2074 m_w[iw].x = xMap[i];
2075 m_w[iw].y = yMap[j];
2076 // Verify the current situation.
2077 if (!WireCheck()) {
2078 std::cerr << m_className << "::ForcesOnWire: Wire " << iw << ".\n"
2079 << " Scan involves a disallowed wire position.\n";
2080 ok = false;
2081 continue;
2082 }
2083 // Recompute the charges for this configuration.
2084 if (!Setup()) {
2085 std::cerr << m_className << "::ForcesOnWire: Wire " << iw << ".\n"
2086 << " Failed to compute charges at a scan point.\n";
2087 ok = false;
2088 continue;
2089 }
2090 // Compute the forces.
2091 double ex = 0., ey = 0.;
2092 ElectricFieldAtWire(iw, ex, ey);
2093 fxMap[i][j] = -ex * wire.e * Internal2Newton;
2094 fyMap[i][j] = -ey * wire.e * Internal2Newton;
2095 }
2096 }
2097 // Place the wire back in its shifted position.
2098 m_w[iw].x = x0;
2099 m_w[iw].y = y0;
2100 Setup();
2101 return ok;
2102}
Garfield::ComponentAnalyticField::ElectricFieldAtWire
bool ElectricFieldAtWire(const unsigned int iw, double &ex, double &ey)
Definition: ComponentAnalyticField.cc:1750
Garfield::ComponentBase::m_debug
bool m_debug
Switch on/off debugging messages.
Definition: ComponentBase.hh:310
Heed::sqrt
DoubleAc sqrt(const DoubleAc &f)
Definition: DoubleAc.cpp:314

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::ComponentBase.

Definition at line 311 of file ComponentAnalyticField.cc.

313 {
314 // If a geometry is present, try to get the bounding box from there.
315 if (m_geometry) {
316 if (m_geometry->GetBoundingBox(x0, y0, z0, x1, y1, z1)) return true;
317 }
318 // Otherwise, return the cell dimensions.
319 if (!m_cellset && !Prepare()) return false;
320 if (m_polar) {
321 double rmax, thetamax;
322 Internal2Polar(m_xmax, m_ymax, rmax, thetamax);
323 x0 = -rmax;
324 y0 = -rmax;
325 x1 = +rmax;
326 y1 = +rmax;
327 } else {
328 x0 = m_xmin;
329 y0 = m_ymin;
330 x1 = m_xmax;
331 y1 = m_ymax;
332 }
333 z0 = m_zmin;
334 z1 = m_zmax;
335 return true;
336}
Garfield::GeometryBase::GetBoundingBox
virtual bool GetBoundingBox(double &xmin, double &ymin, double &zmin, double &xmax, double &ymax, double &zmax)=0
Get the bounding box (envelope of the geometry).

◆ 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 196 of file ComponentAnalyticField.hh.

196 {
197 if (!m_cellset) {
198 if (CellCheck()) CellType();
199 }
200 return GetCellType(m_cellType);
201 }

Referenced by GetCellType(), and PrintCell().

◆ GetGravity()

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

Get the gravity orientation.

Definition at line 1887 of file ComponentAnalyticField.cc.

1888 {
1889
1890 dx = m_down[0];
1891 dy = m_down[1];
1892 dz = m_down[2];
1893}

◆ 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::ComponentBase.

Definition at line 245 of file ComponentAnalyticField.cc.

246 {
247 if (m_geometry) return m_geometry->GetMedium(xin, yin, zin);
248
249 // Make sure the cell is prepared.
250 if (!m_cellset && !Prepare()) return nullptr;
251
252 double xpos = xin, ypos = yin;
253 if (m_polar) Cartesian2Internal(xin, yin, xpos, ypos);
254 // In case of periodicity, move the point into the basic cell.
255 if (m_perx) {
256 xpos -= m_sx * int(round(xin / m_sx));
257 }
258 double arot = 0.;
259 if (m_pery && m_tube) {
260 Cartesian2Polar(xin, yin, xpos, ypos);
261 arot = RadToDegree * m_sy * int(round(DegreeToRad * ypos / m_sy));
262 ypos -= arot;
263 Polar2Cartesian(xpos, ypos, xpos, ypos);
264 } else if (m_pery) {
265 ypos -= m_sy * int(round(ypos / m_sy));
266 }
267
268 // Move the point to the correct side of the plane.
269 if (m_perx && m_ynplan[0] && xpos <= m_coplan[0]) xpos += m_sx;
270 if (m_perx && m_ynplan[1] && xpos >= m_coplan[1]) xpos -= m_sx;
271 if (m_pery && m_ynplan[2] && ypos <= m_coplan[2]) ypos += m_sy;
272 if (m_pery && m_ynplan[3] && ypos >= m_coplan[3]) ypos -= m_sy;
273
274 // In case (xpos, ypos) is located behind a plane there is no field.
275 if (m_tube) {
276 if (!InTube(xpos, ypos, m_cotube, m_ntube)) return nullptr;
277 } else {
278 if ((m_ynplan[0] && xpos < m_coplan[0]) ||
279 (m_ynplan[1] && xpos > m_coplan[1]) ||
280 (m_ynplan[2] && ypos < m_coplan[2]) ||
281 (m_ynplan[3] && ypos > m_coplan[3])) {
282 return nullptr;
283 }
284 }
285
286 // If (xpos, ypos) is within a wire, there is no field either.
287 for (const auto& wire : m_w) {
288 double dx = xpos - wire.x;
289 double dy = ypos - wire.y;
290 // Correct for periodicities.
291 if (m_perx) dx -= m_sx * int(round(dx / m_sx));
292 if (m_pery) dy -= m_sy * int(round(dy / m_sy));
293 // Check the actual position.
294 if (dx * dx + dy * dy < wire.r * wire.r) return nullptr;
295 }
296 return m_medium;
297}

◆ GetNumberOfPlanesPhi()

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

Get the number of equipotential planes at constant phi.

Definition at line 1645 of file ComponentAnalyticField.cc.

1645 {
1646 if (!m_polar) {
1647 return 0;
1648 } else if (m_ynplan[2] && m_ynplan[3]) {
1649 return 2;
1650 } else if (m_ynplan[2] || m_ynplan[3]) {
1651 return 1;
1652 }
1653 return 0;
1654}

◆ GetNumberOfPlanesR()

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

Get the number of equipotential planes at constant radius.

Definition at line 1634 of file ComponentAnalyticField.cc.

1634 {
1635 if (!m_polar) {
1636 return 0;
1637 } else if (m_ynplan[0] && m_ynplan[1]) {
1638 return 2;
1639 } else if (m_ynplan[0] || m_ynplan[1]) {
1640 return 1;
1641 }
1642 return 0;
1643}

◆ GetNumberOfPlanesX()

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

Get the number of equipotential planes at constant x.

Definition at line 1612 of file ComponentAnalyticField.cc.

1612 {
1613 if (m_polar) {
1614 return 0;
1615 } else if (m_ynplan[0] && m_ynplan[1]) {
1616 return 2;
1617 } else if (m_ynplan[0] || m_ynplan[1]) {
1618 return 1;
1619 }
1620 return 0;
1621}

◆ GetNumberOfPlanesY()

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

Get the number of equipotential planes at constant y.

Definition at line 1623 of file ComponentAnalyticField.cc.

1623 {
1624 if (m_polar) {
1625 return 0;
1626 } else if (m_ynplan[2] && m_ynplan[3]) {
1627 return 2;
1628 } else if (m_ynplan[2] || m_ynplan[3]) {
1629 return 1;
1630 }
1631 return 0;
1632}

◆ GetNumberOfWires()

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

Get the number of wires.

Definition at line 226 of file ComponentAnalyticField.hh.

226{ return m_nWires; }

◆ GetPeriodicityPhi()

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

Get the periodicity [degree] in phi.

Definition at line 1493 of file ComponentAnalyticField.cc.

1493 {
1494 if (!m_periodic[1] || (!m_polar && !m_tube)) {
1495 s = 0.;
1496 return false;
1497 }
1498
1499 s = RadToDegree * m_sy;
1500 return true;
1501}
Garfield::ComponentBase::m_periodic
std::array< bool, 3 > m_periodic
Simple periodicity in x, y, z.
Definition: ComponentBase.hh:299

◆ GetPeriodicityX()

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

Get the periodic length in the x-direction.

Definition at line 1473 of file ComponentAnalyticField.cc.

1473 {
1474 if (!m_periodic[0] || m_polar) {
1475 s = 0.;
1476 return false;
1477 }
1478
1479 s = m_sx;
1480 return true;
1481}

◆ GetPeriodicityY()

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

Get the periodic length in the y-direction.

Definition at line 1483 of file ComponentAnalyticField.cc.

1483 {
1484 if (!m_periodic[1] || m_polar) {
1485 s = 0.;
1486 return false;
1487 }
1488
1489 s = m_sy;
1490 return true;
1491}

◆ 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 1726 of file ComponentAnalyticField.cc.

1728 {
1729 if (!m_polar || i >= 2 || (i == 1 && !m_ynplan[3])) {
1730 std::cerr << m_className << "::GetPlanePhi: Index out of range.\n";
1731 return false;
1732 }
1733
1734 phi = RadToDegree * m_coplan[i + 2];
1735 voltage = m_vtplan[i + 2];
1736 label = m_planes[i + 2].type;
1737 return true;
1738}

◆ 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 1712 of file ComponentAnalyticField.cc.

1714 {
1715 if (!m_polar || i >= 2 || (i == 1 && !m_ynplan[1])) {
1716 std::cerr << m_className << "::GetPlaneR: Index out of range.\n";
1717 return false;
1718 }
1719
1720 r = exp(m_coplan[i]);
1721 voltage = m_vtplan[i];
1722 label = m_planes[i].type;
1723 return true;
1724}

◆ 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 1684 of file ComponentAnalyticField.cc.

1686 {
1687 if (m_polar || i >= 2 || (i == 1 && !m_ynplan[1])) {
1688 std::cerr << m_className << "::GetPlaneX: Index out of range.\n";
1689 return false;
1690 }
1691
1692 x = m_coplan[i];
1693 voltage = m_vtplan[i];
1694 label = m_planes[i].type;
1695 return true;
1696}

◆ 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 1698 of file ComponentAnalyticField.cc.

1700 {
1701 if (m_polar || i >= 2 || (i == 1 && !m_ynplan[3])) {
1702 std::cerr << m_className << "::GetPlaneY: Index out of range.\n";
1703 return false;
1704 }
1705
1706 y = m_coplan[i + 2];
1707 voltage = m_vtplan[i + 2];
1708 label = m_planes[i + 2].type;
1709 return true;
1710}

◆ GetTube()

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

Retrieve the tube parameters.

Definition at line 1740 of file ComponentAnalyticField.cc.

1741 {
1742 if (!m_tube) return false;
1743 r = m_cotube;
1744 voltage = m_vttube;
1745 nEdges = m_ntube;
1746 label = m_planes[4].type;
1747 return true;
1748}

◆ GetVoltageRange()

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

Calculate the voltage range [V].

Implements Garfield::ComponentBase.

Definition at line 299 of file ComponentAnalyticField.cc.

299 {
300 // Make sure the cell is prepared.
301 if (!m_cellset && !Prepare()) {
302 std::cerr << m_className << "::GetVoltageRange: Cell not set up.\n";
303 return false;
304 }
305
306 pmin = m_vmin;
307 pmax = m_vmax;
308 return true;
309}

◆ 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 1656 of file ComponentAnalyticField.cc.

1659 {
1660 if (i >= m_nWires) {
1661 std::cerr << m_className << "::GetWire: Index out of range.\n";
1662 return false;
1663 }
1664
1665 if (m_polar) {
1666 double r = 0., theta = 0.;
1667 Internal2Polar(m_w[i].x, m_w[i].y, r, theta);
1668 x = r;
1669 y = theta;
1670 diameter = 2 * m_w[i].r * r;
1671 } else {
1672 x = m_w[i].x;
1673 y = m_w[i].y;
1674 diameter = 2 * m_w[i].r;
1675 }
1676 voltage = m_w[i].v;
1677 label = m_w[i].type;
1678 length = m_w[i].u;
1679 charge = m_w[i].e;
1680 ntrap = m_w[i].nTrap;
1681 return true;
1682}

◆ IsInTrapRadius()

bool Garfield::ComponentAnalyticField::IsInTrapRadius ( 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::ComponentBase.

Definition at line 737 of file ComponentAnalyticField.cc.

740 {
741
742 double x0 = xin;
743 double y0 = yin;
744 if (m_polar) {
745 Cartesian2Internal(xin, yin, x0, y0);
746 }
747 // In case of periodicity, move the point into the basic cell.
748 int nX = 0, nY = 0, nPhi = 0;
749 if (m_perx) {
750 nX = int(round(x0 / m_sx));
751 x0 -= m_sx * nX;
752 }
753 if (m_pery && m_tube) {
754 Cartesian2Polar(xin, yin, x0, y0);
755 nPhi = int(round(DegreeToRad * y0 / m_sy));
756 y0 -= RadToDegree * m_sy * nPhi;
757 Polar2Cartesian(x0, y0, x0, y0);
758 } else if (m_pery) {
759 nY = int(round(y0 / m_sy));
760 y0 -= m_sy * nY;
761 }
762 // Move the point to the correct side of the plane.
763 std::array<bool, 4> shift = {false, false, false, false};
764 if (m_perx && m_ynplan[0] && x0 <= m_coplan[0]) {
765 x0 += m_sx;
766 shift[0] = true;
767 }
768 if (m_perx && m_ynplan[1] && x0 >= m_coplan[1]) {
769 x0 -= m_sx;
770 shift[1] = true;
771 }
772 if (m_pery && m_ynplan[2] && y0 <= m_coplan[2]) {
773 y0 += m_sy;
774 shift[2] = true;
775 }
776 if (m_pery && m_ynplan[3] && y0 >= m_coplan[3]) {
777 y0 -= m_sy;
778 shift[3] = true;
779 }
780
781 for (const auto& wire : m_w) {
782 // Skip wires with the wrong charge.
783 if (qin * wire.e > 0.) continue;
784 const double dxw0 = wire.x - x0;
785 const double dyw0 = wire.y - y0;
786 const double r2 = dxw0 * dxw0 + dyw0 * dyw0;
787 const double rTrap = wire.r * wire.nTrap;
788 if (r2 < rTrap * rTrap) {
789 xw = wire.x;
790 yw = wire.y;
791 rw = wire.r;
792 if (shift[0]) xw -= m_sx;
793 if (shift[1]) xw += m_sx;
794 if (shift[2]) yw -= m_sy;
795 if (shift[3]) yw += m_sy;
796 if (m_pery && m_tube) {
797 double rhow, phiw;
798 Cartesian2Polar(xw, yw, rhow, phiw);
799 phiw += RadToDegree * m_sy * nPhi;
800 Polar2Cartesian(rhow, phiw, xw, yw);
801 } else if (m_pery) {
802 y0 += m_sy * nY;
803 }
804 if (m_perx) xw += m_sx * nX;
805 if (m_polar) {
806 Internal2Cartesian(xw, yw, xw, yw);
807 rw *= exp(wire.x);
808 }
809 if (m_debug) {
810 std::cout << m_className << "::IsInTrapRadius: (" << xin << ", "
811 << yin << ", " << zin << ")" << " within trap radius.\n";
812 }
813 return true;
814 }
815 }
816
817 return false;
818}

◆ IsPolar()

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

Are polar coordinates being used?

Definition at line 164 of file ComponentAnalyticField.hh.

164{ return m_polar; }

◆ IsWireCrossed()

bool Garfield::ComponentAnalyticField::IsWireCrossed ( 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::ComponentBase.

Definition at line 636 of file ComponentAnalyticField.cc.

640 {
641 xc = xx0;
642 yc = yy0;
643 zc = z0;
644
645 if (m_w.empty()) return false;
646
647 double x0 = xx0;
648 double y0 = yy0;
649 double x1 = xx1;
650 double y1 = yy1;
651 if (m_polar) {
652 Cartesian2Internal(xx0, yy0, x0, y0);
653 Cartesian2Internal(xx1, yy1, x1, y1);
654 }
655 const double dx = x1 - x0;
656 const double dy = y1 - y0;
657 const double d2 = dx * dx + dy * dy;
658 // Check that the step length is non-zero.
659 if (d2 < Small) return false;
660 const double invd2 = 1. / d2;
661
662 // Check if a whole period has been crossed.
663 if ((m_perx && fabs(dx) >= m_sx) || (m_pery && fabs(dy) >= m_sy)) {
664 std::cerr << m_className << "::IsWireCrossed:\n"
665 << " Particle crossed more than one period.\n";
666 return false;
667 }
668
669 // Both coordinates are assumed to be located inside
670 // the drift area and inside a drift medium.
671 // This should have been checked before this call.
672
673 const double xm = 0.5 * (x0 + x1);
674 const double ym = 0.5 * (y0 + y1);
675 double dMin2 = 0.;
676 for (const auto& wire : m_w) {
677 double xw = wire.x;
678 if (m_perx) {
679 xw += m_sx * int(round((xm - xw) / m_sx));
680 }
681 double yw = wire.y;
682 if (m_pery) {
683 yw += m_sy * int(round((ym - yw) / m_sy));
684 }
685 // Calculate the smallest distance between track and wire.
686 const double xIn0 = dx * (xw - x0) + dy * (yw - y0);
687 // Check if the minimum is located before (x0, y0).
688 if (xIn0 < 0.) continue;
689 const double xIn1 = -(dx * (xw - x1) + dy * (yw - y1));
690 // Check if the minimum is located behind (x1, y1).
691 if (xIn1 < 0.) continue;
692 // Minimum is located between (x0, y0) and (x1, y1).
693 const double xw0 = xw - x0;
694 const double xw1 = xw - x1;
695 const double yw0 = yw - y0;
696 const double yw1 = yw - y1;
697 const double dw02 = xw0 * xw0 + yw0 * yw0;
698 const double dw12 = xw1 * xw1 + yw1 * yw1;
699 if (xIn1 * xIn1 * dw02 > xIn0 * xIn0 * dw12) {
700 dMin2 = dw02 - xIn0 * xIn0 * invd2;
701 } else {
702 dMin2 = dw12 - xIn1 * xIn1 * invd2;
703 }
704 // Add in the times nTrap to account for the trap radius.
705 const double r2 = wire.r * wire.r;
706 if (dMin2 < r2) {
707 // Wire has been crossed.
708 if (centre) {
709 if (m_polar) {
710 Internal2Cartesian(xw, yw, xc, yc);
711 } else {
712 xc = xw;
713 yc = yw;
714 }
715 } else {
716 // Find the point of intersection.
717 const double p = -xIn0 * invd2;
718 const double q = (dw02 - r2) * invd2;
719 const double s = sqrt(p * p - q);
720 const double t = std::min(-p + s, -p - s);
721 if (m_polar) {
722 Internal2Cartesian(x0 + t * dx, y0 + t * dy, xc, yc);
723 } else {
724 xc = x0 + t * dx;
725 yc = y0 + t * dy;
726 }
727 zc = z0 + t * (z1 - z0);
728 }
729 rc = wire.r;
730 if (m_polar) rc *= exp(wire.x);
731 return true;
732 }
733 }
734 return false;
735}

◆ 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 9783 of file ComponentAnalyticField.cc.

9786 {
9787
9788 //-----------------------------------------------------------------------
9789 // EFMWIR - Computes the dipole moment of a given wire.
9790 //-----------------------------------------------------------------------
9791 if (!m_cellset && !Prepare()) return false;
9792 // Check input parameters.
9793 if (iw >= m_nWires) {
9794 std::cerr << m_className << "::MultipoleMoments:\n"
9795 << " Wire index out of range.\n";
9796 return false;
9797 }
9798 if (eps <= 0.) {
9799 std::cerr << m_className << "::MultipoleMoments:\n"
9800 << " Epsilon must be positive.\n";
9801 return false;
9802 }
9803 if (nPoles == 0) {
9804 std::cerr << m_className << "::MultipoleMoments:\n"
9805 << " Multipole order out of range.\n";
9806 return false;
9807 }
9808 if (rmult <= 0.) {
9809 std::cerr << m_className << "::MultipoleMoments:\n"
9810 << " Radius multiplication factor out of range.\n";
9811 return false;
9812 }
9813
9814 const double xw = m_w[iw].x;
9815 const double yw = m_w[iw].y;
9816 // Set the radius of the wire to 0.
9817 const double rw = m_w[iw].r;
9818 m_w[iw].r = 0.;
9819
9820 // Loop around the wire.
9821 constexpr unsigned int nPoints = 20000;
9822 std::vector<double> angle(nPoints, 0.);
9823 std::vector<double> volt(nPoints, 0.);
9824 std::vector<double> weight(nPoints, 1.);
9825 for (unsigned int i = 0; i < nPoints; ++i) {
9826 // Set angle around wire.
9827 angle[i] = TwoPi * (i + 1.) / nPoints;
9828 // Compute E field, make sure the point is in a free region.
9829 const double x = xw + rmult * rw * cos(angle[i]);
9830 const double y = yw + rmult * rw * sin(angle[i]);
9831 double ex = 0., ey = 0., ez = 0., v = 0.;
9832 if (Field(x, y, 0., ex, ey, ez, v, true) != 0) {
9833 std::cerr << m_className << "::MultipoleMoments:\n"
9834 << " Unexpected location code. Computation stopped.\n";
9835 m_w[iw].r = rw;
9836 return false;
9837 }
9838 // Assign the result to the fitting array.
9839 volt[i] = v;
9840 }
9841 // Restore the wire diameter.
9842 m_w[iw].r = rw;
9843
9844 // Determine the maximum, minimum and average.
9845 double vmin = *std::min_element(volt.cbegin(), volt.cend());
9846 double vmax = *std::max_element(volt.cbegin(), volt.cend());
9847 double vave = std::accumulate(volt.cbegin(), volt.cend(), 0.) / nPoints;
9848 // Subtract the wire potential to centre the data more or less.
9849 for (unsigned int i = 0; i < nPoints; ++i) volt[i] -= vave;
9850 vmax -= vave;
9851 vmin -= vave;
9852
9853 // Perform the fit.
9854 const double vm = 0.5 * fabs(vmin) + fabs(vmax);
9855 double chi2 = 1.e-6 * nPoints * vm * vm;
9856 const double dist = 1.e-3 * (1. + vm);
9857 const unsigned int nPar = 2 * nPoles + 1;
9858 std::vector<double> pars(nPar, 0.);
9859 std::vector<double> epar(nPar, 0.);
9860 pars[0] = 0.5 * (vmax + vmin);
9861 for (unsigned int i = 1; i <= nPoles; ++i) {
9862 pars[2 * i - 1] = 0.5 * (vmax - vmin);
9863 pars[2 * i] = 0.;
9864 }
9865
9866 auto f = [nPoles](const double x, const std::vector<double>& par) {
9867 // EFMFUN
9868 // Sum the series, initial value is the monopole term.
9869 double sum = par[0];
9870 for (unsigned int k = 1; k <= nPoles; ++k) {
9871 // Obtain the Legendre polynomial of this order and add to the series.
9872 const float cphi = cos(x - par[2 * k]);
9873 sum += par[2 * k - 1] * sqrt(k + 0.5) * Numerics::Legendre(k, cphi);
9874 }
9875 return sum;
9876 };
9877
9878 if (!Numerics::LeastSquaresFit(f, pars, epar, angle, volt, weight,
9879 nMaxIter, dist, chi2, eps, m_debug, print)) {
9880 std::cerr << m_className << "::MultipoleMoments:\n"
9881 << " Fitting the multipoles failed; computation stopped.\n";
9882 }
9883 // Plot the result of the fit.
9884 if (plot) {
9885 TCanvas* cfit = new TCanvas();
9886 cfit->SetGridx();
9887 cfit->SetGridy();
9888 cfit->DrawFrame(0., vmin, TwoPi, vmax,
9889 ";Angle around the wire [rad]; Potential - average [V]");
9890 TGraph graph;
9891 graph.SetLineWidth(2);
9892 graph.SetLineColor(kBlack);
9893 graph.DrawGraph(angle.size(), angle.data(), volt.data(), "lsame");
9894 // Sum of contributions.
9895 constexpr unsigned int nP = 1000;
9896 std::array<double, nP> xp;
9897 std::array<double, nP> yp;
9898 for (unsigned int i = 0; i < nP; ++i) {
9899 xp[i] = TwoPi * (i + 1.) / nP;
9900 yp[i] = f(xp[i], pars);
9901 }
9902 graph.SetLineColor(kViolet + 3);
9903 graph.DrawGraph(nP, xp.data(), yp.data(), "lsame");
9904 // Individual contributions.
9905 std::vector<double> parres = pars;
9906 for (unsigned int i = 1; i <= nPoles; ++i) parres[2 * i - 1] = 0.;
9907 for (unsigned int j = 1; j <= nPoles; ++j) {
9908 parres[2 * j - 1] = pars[2 * j - 1];
9909 for (unsigned int i = 0; i < nP; ++i) {
9910 yp[i] = f(xp[i], parres);
9911 }
9912 parres[2 * j - 1] = 0.;
9913 graph.SetLineColor(kAzure + j);
9914 graph.DrawGraph(nP, xp.data(), yp.data(), "lsame");
9915 }
9916 }
9917
9918 // Print the results.
9919 std::cout << m_className << "::MultipoleMoments:\n"
9920 << " Multipole moments for wire " << iw << ":\n"
9921 << " Moment Value Angle [degree]\n";
9922 std::printf(" %6u %15.8f Arbitrary\n", 0, vave);
9923 for (unsigned int i = 1; i <= nPoles; ++i) {
9924 // Remove radial term from the multipole moment.
9925 const double val = pow(rmult * rw, i) * pars[2 * i - 1];
9926 const double phi = RadToDegree * fmod(pars[2 * i], Pi);
9927 std::printf(" %6u %15.8f %15.8f\n", i, val, phi);
9928 }
9929 return true;
9930}
Garfield::Numerics::Legendre
double Legendre(const unsigned int n, const double x)
Legendre polynomials.
Definition: Numerics.hh:119
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:1838
Heed::pow
DoubleAc pow(const DoubleAc &f, double p)
Definition: DoubleAc.cpp:337

◆ PrintCell()

void Garfield::ComponentAnalyticField::PrintCell ( )

Print all available information on the cell.

Definition at line 338 of file ComponentAnalyticField.cc.

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

◆ PrintCharges()

void Garfield::ComponentAnalyticField::PrintCharges ( ) const

Print a list of the point charges.

Definition at line 1598 of file ComponentAnalyticField.cc.

1598 {
1599 std::cout << m_className << "::PrintCharges:\n";
1600 if (m_ch3d.empty()) {
1601 std::cout << " No charges present.\n";
1602 return;
1603 }
1604 std::cout << " x [cm] y [cm] z [cm] charge [fC]\n";
1605 for (const auto& charge : m_ch3d) {
1606 std::cout << " " << std::setw(9) << charge.x << " " << std::setw(9)
1607 << charge.y << " " << std::setw(9) << charge.z << " "
1608 << std::setw(11) << charge.e * FourPiEpsilon0 << "\n";
1609 }
1610}

◆ SetCartesianCoordinates()

void Garfield::ComponentAnalyticField::SetCartesianCoordinates ( )

Use Cartesian coordinates (default).

Definition at line 1560 of file ComponentAnalyticField.cc.

1560 {
1561 if (m_polar) {
1562 std::cout << m_className << "::SetCartesianCoordinates:\n "
1563 << "Switching to Cartesian coordinates; resetting the cell.\n";
1564 CellInit();
1565 }
1566 m_polar = false;
1567}

◆ SetGravity()

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

Set the gravity orientation.

Definition at line 1873 of file ComponentAnalyticField.cc.

1874 {
1875
1876 const double d = sqrt(dx * dx + dy * dy + dz * dz);
1877 if (d > 0.) {
1878 m_down[0] = dx / d;
1879 m_down[1] = dy / d;
1880 m_down[2] = dz / d;
1881 } else {
1882 std::cerr << m_className << "::SetGravity:\n"
1883 << " The gravity vector has zero norm ; ignored.\n";
1884 }
1885}

◆ SetMedium()

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

Set the medium inside the cell.

Definition at line 90 of file ComponentAnalyticField.hh.

90{ m_medium = medium; }

◆ 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 309 of file ComponentAnalyticField.hh.

309{ 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 2104 of file ComponentAnalyticField.cc.

2104 {
2105 if (n == 0) {
2106 std::cerr << m_className << "::SetNumberOfSteps:\n"
2107 << " Number of steps must be > 0.\n";
2108 return;
2109 }
2110 m_nSteps = n;
2111}

◆ SetPeriodicityPhi()

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

Set the periodicity [degree] in phi.

Definition at line 1455 of file ComponentAnalyticField.cc.

1455 {
1456 if (!m_polar && !m_tube) {
1457 std::cerr << m_className << "::SetPeriodicityPhi:\n"
1458 << " Cannot use phi-periodicity with Cartesian coordinates.\n";
1459 return;
1460 }
1461 if (std::abs(360. - s * int(360. / s)) > 1.e-4) {
1462 std::cerr << m_className << "::SetPeriodicityPhi:\n"
1463 << " Phi periods must divide 360; ignored.\n";
1464 return;
1465 }
1466
1467 m_periodic[1] = true;
1468 m_sy = DegreeToRad * s;
1469 m_mtube = int(360. / s);
1470 UpdatePeriodicity();
1471}

◆ SetPeriodicityX()

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

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

Definition at line 1421 of file ComponentAnalyticField.cc.

1421 {
1422 if (m_polar) {
1423 std::cerr << m_className << "::SetPeriodicityX:\n"
1424 << " Cannot use x-periodicity with polar coordinates.\n";
1425 return;
1426 }
1427 if (s < Small) {
1428 std::cerr << m_className << "::SetPeriodicityX:\n"
1429 << " Periodic length must be greater than zero.\n";
1430 return;
1431 }
1432
1433 m_periodic[0] = true;
1434 m_sx = s;
1435 UpdatePeriodicity();
1436}

◆ SetPeriodicityY()

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

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

Definition at line 1438 of file ComponentAnalyticField.cc.

1438 {
1439 if (m_polar) {
1440 std::cerr << m_className << "::SetPeriodicityY:\n"
1441 << " Cannot use y-periodicity with polar coordinates.\n";
1442 return;
1443 }
1444 if (s < Small) {
1445 std::cerr << m_className << "::SetPeriodicityY:\n"
1446 << " Periodic length must be greater than zero.\n";
1447 return;
1448 }
1449
1450 m_periodic[1] = true;
1451 m_sy = s;
1452 UpdatePeriodicity();
1453}

◆ SetPolarCoordinates()

void Garfield::ComponentAnalyticField::SetPolarCoordinates ( )

Use polar coordinates.

Definition at line 1569 of file ComponentAnalyticField.cc.

1569 {
1570 if (!m_polar) {
1571 std::cout << m_className << "::SetPolarCoordinates:\n "
1572 << "Switching to polar coordinates; resetting the cell.\n";
1573 CellInit();
1574 }
1575 m_polar = true;
1576 // Set default phi period.
1577 m_pery = true;
1578 m_sy = TwoPi;
1579}

◆ 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 1847 of file ComponentAnalyticField.cc.

1850 {
1851 if (std::abs(xmax - xmin) < Small || std::abs(ymax - ymin) < Small) {
1852 std::cerr << m_className << "::SetScanningArea:\n"
1853 << " Zero range not permitted.\n";
1854 return;
1855 }
1856 m_scanRange = ScanningRange::User;
1857 m_xScanMin = std::min(xmin, xmax);
1858 m_xScanMax = std::max(xmin, xmax);
1859 m_yScanMin = std::min(ymin, ymax);
1860 m_yScanMax = std::max(ymin, ymax);
1861}

◆ 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 1863 of file ComponentAnalyticField.cc.

1863 {
1864 m_scanRange = ScanningRange::FirstOrder;
1865 if (scale > 0.) {
1866 m_scaleRange = scale;
1867 } else {
1868 std::cerr << m_className << "::SetScanningAreaFirstOrder:\n"
1869 << " Scaling factor must be > 0.\n";
1870 }
1871}

◆ 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 269 of file ComponentAnalyticField.hh.

269{ 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 1831 of file ComponentAnalyticField.cc.

1832 {
1833 if (nX < 2) {
1834 std::cerr << m_className << "::SetScanningGrid:\n"
1835 << " Number of x-lines must be > 1.\n";
1836 } else {
1837 m_nScanX = nX;
1838 }
1839 if (nY < 2) {
1840 std::cerr << m_className << "::SetScanningGrid:\n"
1841 << " Number of y-lines must be > 1.\n";
1842 } else {
1843 m_nScanY = nY;
1844 }
1845}

◆ 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::ComponentBase.

Definition at line 59 of file ComponentAnalyticField.hh.

61 {
62 wx = wy = wz = 0.;
63 double volt = 0.;
64 if (!m_sigset) PrepareSignals();
65 Wfield(x, y, z, wx, wy, wz, volt, label, false);
66 }

◆ 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::ComponentBase.

Definition at line 67 of file ComponentAnalyticField.hh.

68 {
69 double wx = 0., wy = 0., wz = 0.;
70 double volt = 0.;
71 if (!m_sigset) PrepareSignals();
72 Wfield(x, y, z, wx, wy, wz, volt, label, true);
73 return volt;
74 }

◆ 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 2113 of file ComponentAnalyticField.cc.

2116 {
2117 if (!m_cellset && !Prepare()) {
2118 std::cerr << m_className << "::WireDisplacement: Cell not set up.\n";
2119 return false;
2120 }
2121 if (iw >= m_nWires) {
2122 std::cerr << m_className
2123 << "::WireDisplacement: Wire index out of range.\n";
2124 return false;
2125 }
2126 if (m_polar) {
2127 std::cerr << m_className
2128 << "::WireDisplacement: Cannot handle polar cells.\n";
2129 return false;
2130 }
2131 const auto& wire = m_w[iw];
2132 // Save the original coordinates.
2133 const double x0 = wire.x;
2134 const double y0 = wire.y;
2135 // First-order approximation.
2136 if (!Setup()) {
2137 std::cerr << m_className << "::WireDisplacement:\n"
2138 << " Charge calculation failed at central position.\n";
2139 return false;
2140 }
2141
2142 double fx = 0., fy = 0.;
2143 if (m_useElectrostaticForce) {
2144 double ex = 0., ey = 0.;
2145 ElectricFieldAtWire(iw, ex, ey);
2146 fx -= ex * wire.e * Internal2Newton;
2147 fy -= ey * wire.e * Internal2Newton;
2148 }
2149 if (m_useGravitationalForce) {
2150 // Mass per unit length [kg / cm].
2151 const double m = 1.e-3 * wire.density * Pi * wire.r * wire.r;
2152 fx -= m_down[0] * GravitationalAcceleration * m;
2153 fy -= m_down[1] * GravitationalAcceleration * m;
2154 }
2155 const double u2 = wire.u * wire.u;
2156 const double shiftx =
2157 -125 * fx * u2 / (GravitationalAcceleration * wire.tension);
2158 const double shifty =
2159 -125 * fy * u2 / (GravitationalAcceleration * wire.tension);
2160 // Get the elongation from this.
2161 const double s = 4 * sqrt(shiftx * shiftx + shifty * shifty) / wire.u;
2162 double length = wire.u;
2163 if (s > Small) {
2164 const double t = sqrt(1 + s * s);
2165 length *= 0.5 * (t + log(s + t) / s);
2166 }
2167 stretch = (length - wire.u) / wire.u;
2168 if (print) {
2169 std::cout << m_className
2170 << "::WireDisplacement:\n"
2171 << " Forces and displacement in 0th order.\n"
2172 << " Wire information: number = " << iw << "\n"
2173 << " type = " << wire.type << "\n"
2174 << " location = (" << wire.x << ", " << wire.y
2175 << ") cm\n"
2176 << " voltage = " << wire.v << " V\n"
2177 << " length = " << wire.u << " cm\n"
2178 << " tension = " << wire.tension << " g\n"
2179 << " In this position: Fx = " << fx << " N/cm\n"
2180 << " Fy = " << fy << " N/cm\n"
2181 << " x-shift = " << shiftx << " cm\n"
2182 << " y-shift = " << shifty << " cm\n"
2183 << " stretch = " << 100. * stretch << "%\n";
2184 }
2185 if (!detailed) {
2186 csag = {0.};
2187 xsag = {shiftx};
2188 ysag = {shifty};
2189 return true;
2190 }
2191 std::vector<double> xMap(m_nScanX, 0.);
2192 std::vector<double> yMap(m_nScanY, 0.);
2193 std::vector<std::vector<double> > fxMap(m_nScanX,
2194 std::vector<double>(m_nScanY, 0.));
2195 std::vector<std::vector<double> > fyMap(m_nScanX,
2196 std::vector<double>(m_nScanY, 0.));
2197 if (!ForcesOnWire(iw, xMap, yMap, fxMap, fyMap)) {
2198 std::cerr << m_className << "::WireDisplacement:\n"
2199 << " Could not compute the electrostatic force table.\n";
2200 return false;
2201 }
2202 // Compute the detailed wire shift.
2203 if (!SagDetailed(wire, xMap, yMap, fxMap, fyMap, csag, xsag, ysag)) {
2204 std::cerr << m_className << "::WireDisplacement: Wire " << iw << ".\n"
2205 << " Computation of the wire sag failed.\n";
2206 return false;
2207 }
2208 // Verify that the wire is in range.
2209 const double sxmin = xMap.front();
2210 const double sxmax = xMap.back();
2211 const double symin = yMap.front();
2212 const double symax = yMap.back();
2213 const unsigned int nSag = xsag.size();
2214 bool outside = false;
2215 length = 0.;
2216 double xAvg = 0.;
2217 double yAvg = 0.;
2218 double xMax = 0.;
2219 double yMax = 0.;
2220 for (unsigned int i = 0; i < nSag; ++i) {
2221 if (x0 + xsag[i] < sxmin || x0 + xsag[i] > sxmax || y0 + ysag[i] < symin ||
2222 y0 + ysag[i] > symax) {
2223 outside = true;
2224 }
2225 xAvg += xsag[i];
2226 yAvg += ysag[i];
2227 xMax = std::max(xMax, std::abs(xsag[i]));
2228 yMax = std::max(yMax, std::abs(ysag[i]));
2229 if (i == 0) continue;
2230 const double dx = xsag[i] - xsag[i - 1];
2231 const double dy = ysag[i] - ysag[i - 1];
2232 const double dz = csag[i] - csag[i - 1];
2233 length += sqrt(dx * dx + dy * dy + dz * dz);
2234 }
2235 xAvg /= nSag;
2236 yAvg /= nSag;
2237 stretch = (length - wire.u) / wire.u;
2238 // Warn if a point outside the scanning area was found.
2239 if (outside) {
2240 std::cerr
2241 << m_className << "::WireDisplacement: Warning.\n "
2242 << "The wire profile is located partially outside the scanning area.\n";
2243 }
2244 if (print) {
2245 std::cout << " Sag profile for wire " << iw << ".\n"
2246 << " Point z [cm] x-sag [um] y-sag [um]\n";
2247 for (unsigned int i = 0; i < nSag; ++i) {
2248 std::printf(" %3d %10.4f %10.4f %10.4f\n",
2249 i, csag[i], xsag[i] * 1.e4, ysag[i] * 1.e4);
2250 }
2251 std::printf(" Average sag in x and y: %10.4f and %10.4f micron\n",
2252 1.e4 * xAvg, 1.e4 * yAvg);
2253 std::printf(" Maximum sag in x and y: %10.4f and %10.4f micron\n",
2254 1.e4 * xMax, 1.e4 * yMax);
2255 std::cout << " Elongation: " << 100. * stretch << "%\n";
2256 }
2257 return true;
2258}
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:1895
The documentation for this class was generated from the following files: