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	GetElementaryCell (double &x0, double &y0, double &z0, double &x1, double &y1, double &z1) override
	Get the coordinates of the elementary cell.
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.)
void	AddStripOnPlaneY (const char direction, const double y, const double smin, const double smax, const std::string &label, const double gap=-1.)
	Add a strip in the x or z direction on an existing plane at constant y.
void	AddStripOnPlaneR (const char direction, const double r, const double smin, const double smax, const std::string &label, const double gap=-1.)
	Add a strip in the phi or z direction on an existing plane at constant radius.
void	AddStripOnPlanePhi (const char direction, const double phi, const double smin, const double smax, const std::string &label, const double gap=-1.)
	Add a strip in the r or z direction on an existing plane at constant phi.
void	AddPixelOnPlaneX (const double x, const double ymin, const double ymax, const double zmin, const double zmax, const std::string &label, const double gap=-1., const double rot=0.)
void	AddPixelOnPlaneY (const double y, const double xmin, const double xmax, const double zmin, const double zmax, const std::string &label, const double gap=-1., const double rot=0.)
	Add a pixel on an existing plane at constant y.
void	AddPixelOnPlaneR (const double r, const double phimin, const double phimax, const double zmin, const double zmax, const std::string &label, const double gap=-1.)
	Add a pixel on an existing plane at constant radius.
void	AddPixelOnPlanePhi (const double phi, const double rmin, const double rmax, const double zmin, const double zmax, const std::string &label, const double gap=-1.)
	Add a pixel on an existing plane at constant phi.
void	SetPeriodicityX (const double s)
	Set the periodic length [cm] in the x-direction.
void	SetPeriodicityY (const double s)
	Set the periodic length [cm] in the y-direction.
void	SetPeriodicityPhi (const double phi)
	Set the periodicity [degree] in phi.
bool	GetPeriodicityX (double &s)
	Get the periodic length in the x-direction.
bool	GetPeriodicityY (double &s)
	Get the periodic length in the y-direction.
bool	GetPeriodicityPhi (double &s)
	Get the periodicity [degree] in phi.
void	SetCartesianCoordinates ()
	Use Cartesian coordinates (default).
void	SetPolarCoordinates ()
	Use polar coordinates.
bool	IsPolar () const
	Are polar coordinates being used?
void	PrintCell ()
	Print all available information on the cell.
void	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::Component
	Component ()=delete
	Default constructor.
	Component (const std::string &name)
	Constructor.
virtual	~Component ()
	Destructor.
virtual void	SetGeometry (Geometry *geo)
	Define the geometry.
virtual void	Clear ()
	Reset.
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 double	DelayedWeightingPotential (const double x, const double y, const double z, const double t, const std::string &label)
virtual void	MagneticField (const double x, const double y, const double z, double &bx, double &by, double &bz, int &status)
void	SetMagneticField (const double bx, const double by, const double bz)
	Set a constant magnetic field.
virtual bool	IsReady ()
	Ready for use?
virtual bool	GetBoundingBox (double &xmin, double &ymin, double &zmin, double &xmax, double &ymax, double &zmax)
	Get the bounding box coordinates.
virtual bool	GetElementaryCell (double &xmin, double &ymin, double &zmin, double &xmax, double &ymax, double &zmax)
	Get the coordinates of the elementary cell.
double	IntegrateFluxCircle (const double xc, const double yc, const double r, const unsigned int nI=50)
double	IntegrateFluxSphere (const double xc, const double yc, const double zc, const double r, const unsigned int nI=20)
double	IntegrateFluxParallelogram (const double x0, const double y0, const double z0, const double dx1, const double dy1, const double dz1, const double dx2, const double dy2, const double dz2, const unsigned int nU=20, const unsigned int nV=20)
double	IntegrateWeightingFluxParallelogram (const std::string &label, const double x0, const double y0, const double z0, const double dx1, const double dy1, const double dz1, const double dx2, const double dy2, const double dz2, const unsigned int nU=20, const unsigned int nV=20)
double	IntegrateFluxLine (const double x0, const double y0, const double z0, const double x1, const double y1, const double z1, const double xp, const double yp, const double zp, const unsigned int nI, const int isign=0)
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 direction.
void	EnablePeriodicityY (const bool on=true)
	Enable simple periodicity in the direction.
void	EnablePeriodicityZ (const bool on=true)
	Enable simple periodicity in the direction.
void	IsPeriodic (bool &perx, bool &pery, bool &perz)
	Return periodicity flags.
void	EnableMirrorPeriodicityX (const bool on=true)
	Enable mirror periodicity in the direction.
void	EnableMirrorPeriodicityY (const bool on=true)
	Enable mirror periodicity in the direction.
void	EnableMirrorPeriodicityZ (const bool on=true)
	Enable mirror periodicity in the direction.
void	IsMirrorPeriodic (bool &perx, bool &pery, bool &perz)
	Return mirror periodicity flags.
void	EnableAxialPeriodicityX (const bool on=true)
	Enable axial periodicity in the direction.
void	EnableAxialPeriodicityY (const bool on=true)
	Enable axial periodicity in the direction.
void	EnableAxialPeriodicityZ (const bool on=true)
	Enable axial periodicity in the direction.
void	IsAxiallyPeriodic (bool &perx, bool &pery, bool &perz)
	Return axial periodicity flags.
void	EnableRotationSymmetryX (const bool on=true)
	Enable rotation symmetry around the axis.
void	EnableRotationSymmetryY (const bool on=true)
	Enable rotation symmetry around the axis.
void	EnableRotationSymmetryZ (const bool on=true)
	Enable rotation symmetry around the axis.
void	IsRotationSymmetric (bool &rotx, bool &roty, bool &rotz)
	Return rotation symmetry flags.
void	EnableDebugging ()
	Switch on debugging messages.
void	DisableDebugging ()
	Switch off debugging messages.
virtual bool	HasAttachmentMap () const
	Does the component have attachment maps?
virtual bool	HasVelocityMap () const
	Does the component have velocity maps?
virtual bool	ElectronAttachment (const double, const double, const double, double &eta)
	Get the electron attachment coefficient.
virtual bool	HoleAttachment (const double, const double, const double, double &eta)
	Get the hole attachment coefficient.
virtual bool	ElectronVelocity (const double, const double, const double, double &vx, double &vy, double &vz)
	Get the electron drift velocity.
virtual bool	HoleVelocity (const double, const double, const double, double &vx, double &vy, double &vz)
	Get the hole drift velocity.
virtual bool	GetElectronLifetime (const double, const double, const double, double &etau)
virtual bool	GetHoleLifetime (const double, const double, const double, double &htau)

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

Detailed Description

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

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

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

Constructor & Destructor Documentation

ComponentAnalyticField()

Garfield::ComponentAnalyticField::ComponentAnalyticField

(

)

Constructor.

Definition at line 249 of file ComponentAnalyticField.cc.

                                               : Component("AnalyticField") {
  CellInit();
}

~ComponentAnalyticField()

Garfield::ComponentAnalyticField::~ComponentAnalyticField ( )

inline

Destructor.

Definition at line 21 of file ComponentAnalyticField.hh.

{}

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

                                                                       {
  // Convert from fC to internal units (division by 4 pi epsilon0).
  Charge3d charge;
  charge.x = x;
  charge.y = y;
  charge.z = z;
  charge.e = q / FourPiEpsilon0;
  m_ch3d.push_back(std::move(charge));
}

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

                      {
  if (!m_polar || (!m_ynplan[2] && !m_ynplan[3])) {
    std::cerr << m_className << "::AddPixelOnPlanePhi:\n"
              << "    There are no planes at constant phi.\n";
    return;
  }
 
  if (fabs(rmax - rmin) < Small || fabs(zmax - zmin) < Small) {
    std::cerr << m_className << "::AddPixelOnPlanePhi:\n"
              << "    Pixel width must be greater than zero.\n";
    return;
  }
  if (rmin < Small || rmax < Small) {
    std::cerr << m_className << "::AddPixelOnPlanePhi:\n"
              << "    Radius must be greater than zero.\n";
    return;
  }
  Pixel pixel;
  pixel.type = label;
  pixel.ind = -1;
  const double smin = log(rmin);
  const double smax = log(rmax);
  pixel.smin = std::min(smin, smax);
  pixel.smax = std::max(smin, smax);
  pixel.zmin = std::min(zmin, zmax);
  pixel.zmax = std::max(zmin, zmax);
  pixel.gap = gap > Small ? DegreeToRad * gap : -1.;
 
  int iplane = 2;
  if (m_ynplan[3]) {
    const double d0 = fabs(m_coplan[2] - phi * DegreeToRad);
    const double d1 = fabs(m_coplan[3] - phi * DegreeToRad);
    if (d1 < d0) iplane = 3;
  }
 
  m_planes[iplane].pixels.push_back(std::move(pixel));
}

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

                      {
  if (!m_polar || (!m_ynplan[0] && !m_ynplan[1])) {
    std::cerr << m_className << "::AddPixelOnPlaneR:\n"
              << "    There are no planes at constant r.\n";
    return;
  }
 
  if (fabs(phimax - phimin) < Small || fabs(zmax - zmin) < Small) {
    std::cerr << m_className << "::AddPixelOnPlaneR:\n"
              << "    Pixel width must be greater than zero.\n";
    return;
  }
 
  Pixel pixel;
  pixel.type = label;
  pixel.ind = -1;
  const double smin = phimin * DegreeToRad;
  const double smax = phimax * DegreeToRad;
  pixel.smin = std::min(smin, smax);
  pixel.smax = std::max(smin, smax);
  pixel.zmin = std::min(zmin, zmax);
  pixel.zmax = std::max(zmin, zmax);
  pixel.gap = gap > Small ? gap : -1.;
 
  int iplane = 0;
  if (m_ynplan[1]) {
    const double rho = r > 0. ? log(r) : -25.;
    const double d0 = fabs(m_coplan[0] - rho);
    const double d1 = fabs(m_coplan[1] - rho);
    if (d1 < d0) iplane = 1;
  }
 
  m_planes[iplane].pixels.push_back(std::move(pixel));
}

AddPixelOnPlaneX()

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

Add a pixel on an existing plane at constant x.

Parameters

x	coordinate of the plane.
ymin	lower limit of the pixel cell in y,
ymax	upper limit of the pixel cell in y.
zmin	lower limit of the pixel cell in z.
zmax	upper limit of the pixel cell in z.
label	weighting field identifier.
gap	distance to the opposite plane (optional).
rot	rotation angle (rad) of the pixel (optional).

Definition at line 1281 of file ComponentAnalyticField.cc.

                      {
  if (m_polar || (!m_ynplan[0] && !m_ynplan[1])) {
    std::cerr << m_className << "::AddPixelOnPlaneX:\n"
              << "    There are no planes at constant x.\n";
    return;
  }
 
  if (fabs(ymax - ymin) < Small || fabs(zmax - zmin) < Small) {
    std::cerr << m_className << "::AddPixelOnPlaneX:\n"
              << "    Pixel width must be greater than zero.\n";
    return;
  }
 
  Pixel pixel;
  pixel.type = label;
  pixel.ind = -1;
  pixel.smin = std::min(ymin, ymax);
  pixel.smax = std::max(ymin, ymax);
  pixel.zmin = std::min(zmin, zmax);
  pixel.zmax = std::max(zmin, zmax);
  pixel.gap = gap > Small ? gap : -1.;
  if (fabs(rot) > 1.e-9) {
    pixel.cphi = cos(rot);
    pixel.sphi = sin(rot);
  }
 
  int iplane = 0;
  if (m_ynplan[1]) {
    const double d0 = fabs(m_coplan[0] - x);
    const double d1 = fabs(m_coplan[1] - x);
    if (d1 < d0) iplane = 1;
  }
 
  m_planes[iplane].pixels.push_back(std::move(pixel));
}

AddPixelOnPlaneY()

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

Add a pixel on an existing plane at constant y.

Definition at line 1320 of file ComponentAnalyticField.cc.

                      {
  if (m_polar || (!m_ynplan[2] && !m_ynplan[3])) {
    std::cerr << m_className << "::AddPixelOnPlaneY:\n"
              << "    There are no planes at constant y.\n";
    return;
  }
 
  if (fabs(xmax - xmin) < Small || fabs(zmax - zmin) < Small) {
    std::cerr << m_className << "::AddPixelOnPlaneY:\n"
              << "    Pixel width must be greater than zero.\n";
    return;
  }
 
  Pixel pixel;
  pixel.type = label;
  pixel.ind = -1;
  pixel.smin = std::min(xmin, xmax);
  pixel.smax = std::max(xmin, xmax);
  pixel.zmin = std::min(zmin, zmax);
  pixel.zmax = std::max(zmin, zmax);
  pixel.gap = gap > Small ? gap : -1.;
  if (fabs(rot) > 1.e-9) {
    pixel.cphi = cos(rot);
    pixel.sphi = sin(rot);
  }
 
  int iplane = 2;
  if (m_ynplan[3]) {
    const double d0 = fabs(m_coplan[2] - y);
    const double d1 = fabs(m_coplan[3] - y);
    if (d1 < d0) iplane = 3;
  }
 
  m_planes[iplane].pixels.push_back(std::move(pixel));
}

AddPlanePhi()

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

Add a plane at constant phi.

Definition at line 1043 of file ComponentAnalyticField.cc.

                                                                 {
  if (!m_polar) {
    std::cerr << m_className << "::AddPlanePhi:\n"
              << "    Not compatible with Cartesian coordinates; ignored.\n";
    return;
  }
  if (m_ynplan[2] && m_ynplan[3]) {
    std::cerr << m_className << "::AddPlanePhi:\n"
              << "    Cannot have more than two planes at constant phi.\n";
    return;
  }
 
  if (m_ynplan[2]) {
    m_ynplan[3] = true;
    m_coplan[3] = phi * DegreeToRad;
    m_vtplan[3] = v;
    m_planes[3].type = label;
    m_planes[3].ind = -1;
  } else {
    m_ynplan[2] = true;
    m_coplan[2] = phi * DegreeToRad;
    m_vtplan[2] = v;
    m_planes[2].type = label;
    m_planes[2].ind = -1;
    // Switch off default periodicity.
    if (m_pery && std::abs(m_sy - TwoPi) < 1.e-4) {
      m_pery = false;
    }
  }
 
  // Force recalculation of the capacitance and signal matrices.
  m_cellset = false;
  m_sigset = false;
}

AddPlaneR()

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

Add a plane at constant radius.

Definition at line 1005 of file ComponentAnalyticField.cc.

                                                               {
  if (!m_polar) {
    std::cerr << m_className << "::AddPlaneR:\n"
              << "    Not compatible with Cartesian coordinates; ignored.\n";
    return;
  }
  if (r <= 0.) {
    std::cerr << m_className << "::AddPlaneR:\n"
              << "    Radius must be larger than zero; plane ignored.\n";
    return;
  }
 
  if (m_ynplan[0] && m_ynplan[1]) {
    std::cerr << m_className << "::AddPlaneR:\n"
              << "    Cannot have more than two circular planes.\n";
    return;
  }
 
  if (m_ynplan[0]) {
    m_ynplan[1] = true;
    m_coplan[1] = log(r);
    m_vtplan[1] = v;
    m_planes[1].type = label;
    m_planes[1].ind = -1;
  } else {
    m_ynplan[0] = true;
    m_coplan[0] = log(r);
    m_vtplan[0] = v;
    m_planes[0].type = label;
    m_planes[0].ind = -1;
  }
 
  // Force recalculation of the capacitance and signal matrices.
  m_cellset = false;
  m_sigset = false;
}

AddPlaneX()

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

Add a plane at constant x.

Definition at line 941 of file ComponentAnalyticField.cc.

                                                               {
  if (m_polar) {
    std::cerr << m_className << "::AddPlaneX:\n"
              << "    Not compatible with polar coordinates; ignored.\n";
    return;
  }
  if (m_ynplan[0] && m_ynplan[1]) {
    std::cerr << m_className << "::AddPlaneX:\n"
              << "    Cannot have more than two planes at constant x.\n";
    return;
  }
 
  if (m_ynplan[0]) {
    m_ynplan[1] = true;
    m_coplan[1] = x;
    m_vtplan[1] = v;
    m_planes[1].type = label;
    m_planes[1].ind = -1;
  } else {
    m_ynplan[0] = true;
    m_coplan[0] = x;
    m_vtplan[0] = v;
    m_planes[0].type = label;
    m_planes[0].ind = -1;
  }
 
  // Force recalculation of the capacitance and signal matrices.
  m_cellset = false;
  m_sigset = false;
}

AddPlaneY()

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

Add a plane at constant y.

Definition at line 973 of file ComponentAnalyticField.cc.

                                                               {
  if (m_polar) {
    std::cerr << m_className << "::AddPlaneY:\n"
              << "    Not compatible with polar coordinates; ignored.\n";
    return;
  }
  if (m_ynplan[2] && m_ynplan[3]) {
    std::cerr << m_className << "::AddPlaneY:\n"
              << "    Cannot have more than two planes at constant y.\n";
    return;
  }
 
  if (m_ynplan[2]) {
    m_ynplan[3] = true;
    m_coplan[3] = y;
    m_vtplan[3] = v;
    m_planes[3].type = label;
    m_planes[3].ind = -1;
  } else {
    m_ynplan[2] = true;
    m_coplan[2] = y;
    m_vtplan[2] = v;
    m_planes[2].type = label;
    m_planes[2].ind = -1;
  }
 
  // Force recalculation of the capacitance and signal matrices.
  m_cellset = false;
  m_sigset = false;
}

Referenced by main().

AddReadout()

void Garfield::ComponentAnalyticField::AddReadout

(

const std::string &

label

)

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

Definition at line 3535 of file ComponentAnalyticField.cc.

                                                              {
  // Check if this readout group already exists.
  if (std::find(m_readout.begin(), m_readout.end(), label) != m_readout.end()) {
    std::cout << m_className << "::AddReadout:\n";
    std::cout << "    Readout group " << label << " already exists.\n";
    return;
  }
  m_readout.push_back(label);
 
  unsigned int nWiresFound = 0;
  for (const auto& wire : m_w) {
    if (wire.type == label) ++nWiresFound;
  }
 
  unsigned int nPlanesFound = 0;
  unsigned int nStripsFound = 0;
  unsigned int nPixelsFound = 0;
  for (int i = 0; i < 5; ++i) {
    if (m_planes[i].type == label) ++nPlanesFound;
    for (const auto& strip : m_planes[i].strips1) {
      if (strip.type == label) ++nStripsFound;
    }
    for (const auto& strip : m_planes[i].strips2) {
      if (strip.type == label) ++nStripsFound;
    }
    for (const auto& pixel : m_planes[i].pixels) {
      if (pixel.type == label) ++nPixelsFound;
    }
  }
 
  if (nWiresFound == 0 && nPlanesFound == 0 && nStripsFound == 0 &&
      nPixelsFound == 0) {
    std::cerr << m_className << "::AddReadout:\n";
    std::cerr << "    At present there are no wires, planes or strips\n";
    std::cerr << "    associated to readout group " << label << ".\n";
  } else {
    std::cout << m_className << "::AddReadout:\n";
    std::cout << "    Readout group " << label << " comprises:\n";
    if (nWiresFound > 1) {
      std::cout << "      " << nWiresFound << " wires\n";
    } else if (nWiresFound == 1) {
      std::cout << "      1 wire\n";
    }
    if (nPlanesFound > 1) {
      std::cout << "      " << nPlanesFound << " planes\n";
    } else if (nPlanesFound == 1) {
      std::cout << "      1 plane\n";
    }
    if (nStripsFound > 1) {
      std::cout << "      " << nStripsFound << " strips\n";
    } else if (nStripsFound == 1) {
      std::cout << "      1 strip\n";
    }
    if (nPixelsFound > 1) {
      std::cout << "      " << nPixelsFound << " pixels\n";
    } else if (nPixelsFound == 1) {
      std::cout << "      1 pixel\n";
    }
  }
 
  m_sigset = false;
}

Referenced by main().

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

                                                                  {
  if (!m_polar || (!m_ynplan[2] && !m_ynplan[3])) {
    std::cerr << m_className << "::AddStripOnPlanePhi:\n"
              << "    There are no planes at constant phi.\n";
    return;
  }
 
  if (dir != 'r' && dir != 'R' && dir != 'z' && dir != 'Z') {
    std::cerr << m_className << "::AddStripOnPlanePhi:\n"
              << "    Invalid direction (" << dir << ").\n"
              << "    Only strips in r or z direction are possible.\n";
    return;
  }
 
  if (fabs(smax - smin) < Small) {
    std::cerr << m_className << "::AddStripOnPlanePhi:\n"
              << "    Strip width must be greater than zero.\n";
    return;
  }
 
  Strip newStrip;
  newStrip.type = label;
  newStrip.ind = -1;
  if (dir== 'z' || dir == 'Z') {
    if (smin < Small || smax < Small) {
      std::cerr << m_className << "::AddStripOnPlanePhi:\n"
                << "    Radius must be greater than zero.\n";
      return;
    }
    const double rhomin = log(smin);
    const double rhomax = log(smax);
    newStrip.smin = std::min(rhomin, rhomax);
    newStrip.smax = std::max(rhomin, rhomax);
  } else {
    newStrip.smin = std::min(smin, smax);
    newStrip.smax = std::max(smin, smax);
  }
  newStrip.gap = gap > Small ? DegreeToRad * gap : -1.;
 
  int iplane = 2;
  if (m_ynplan[3]) {
    const double d2 = fabs(m_coplan[2] - phi * DegreeToRad);
    const double d3 = fabs(m_coplan[3] - phi * DegreeToRad);
    if (d3 < d2) iplane = 3;
  }
 
  if (dir == 'r' || dir == 'R') {
    m_planes[iplane].strips1.push_back(std::move(newStrip));
  } else {
    m_planes[iplane].strips2.push_back(std::move(newStrip));
  }
}

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

                                                                {
  if (!m_polar || (!m_ynplan[0] && !m_ynplan[1])) {
    std::cerr << m_className << "::AddStripOnPlaneR:\n"
              << "    There are no planes at constant r.\n";
    return;
  }
 
  if (dir != 'p' && dir != 'P' && dir != 'z' && dir != 'Z') {
    std::cerr << m_className << "::AddStripOnPlaneR:\n"
              << "    Invalid direction (" << dir << ").\n"
              << "    Only strips in p(hi) or z direction are possible.\n";
    return;
  }
 
  if (fabs(smax - smin) < Small) {
    std::cerr << m_className << "::AddStripOnPlaneR:\n"
              << "    Strip width must be greater than zero.\n";
    return;
  }
 
  Strip newStrip;
  newStrip.type = label;
  newStrip.ind = -1;
  if (dir == 'z' || dir == 'Z') {
    const double phimin = smin * DegreeToRad;
    const double phimax = smax * DegreeToRad;
    newStrip.smin = std::min(phimin, phimax);
    newStrip.smax = std::max(phimin, phimax);
  } else {
    newStrip.smin = std::min(smin, smax);
    newStrip.smax = std::max(smin, smax);
  }
  newStrip.gap = gap > Small ? gap : -1.;
 
  int iplane = 0;
  if (m_ynplan[1]) {
    const double rho = r > 0. ? log(r) : -25.;
    const double d0 = fabs(m_coplan[0] - rho);
    const double d1 = fabs(m_coplan[1] - rho);
    if (d1 < d0) iplane = 1;
  }
 
  if (dir == 'p' || dir == 'P') {
    m_planes[iplane].strips1.push_back(std::move(newStrip));
  } else {
    m_planes[iplane].strips2.push_back(std::move(newStrip));
  }
}

AddStripOnPlaneX()

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

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

Parameters

direction	'y' or 'z'.
x	coordinate of the plane.
smin	lower limit of the strip in y or z.
smax	upper limit of the strip in y or z.
label	weighting field identifier.
gap	distance to the opposite plane (optional).

Definition at line 1079 of file ComponentAnalyticField.cc.

                                                                {
  if (m_polar || (!m_ynplan[0] && !m_ynplan[1])) {
    std::cerr << m_className << "::AddStripOnPlaneX:\n"
              << "    There are no planes at constant x.\n";
    return;
  }
 
  if (dir != 'y' && dir != 'Y' && dir != 'z' && dir != 'Z') {
    std::cerr << m_className << "::AddStripOnPlaneX:\n"
              << "    Invalid direction (" << dir << ").\n"
              << "    Only strips in y or z direction are possible.\n";
    return;
  }
 
  if (fabs(smax - smin) < Small) {
    std::cerr << m_className << "::AddStripOnPlaneX:\n"
              << "    Strip width must be greater than zero.\n";
    return;
  }
 
  Strip newStrip;
  newStrip.type = label;
  newStrip.ind = -1;
  newStrip.smin = std::min(smin, smax);
  newStrip.smax = std::max(smin, smax);
  newStrip.gap = gap > Small ? gap : -1.;
 
  int iplane = 0;
  if (m_ynplan[1]) {
    const double d0 = fabs(m_coplan[0] - x);
    const double d1 = fabs(m_coplan[1] - x);
    if (d1 < d0) iplane = 1;
  }
 
  if (dir == 'y' || dir == 'Y') {
    m_planes[iplane].strips1.push_back(std::move(newStrip));
  } else {
    m_planes[iplane].strips2.push_back(std::move(newStrip));
  }
}

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

                                                                {
  if (m_polar || (!m_ynplan[2] && !m_ynplan[3])) {
    std::cerr << m_className << "::AddStripOnPlaneY:\n"
              << "    There are no planes at constant y.\n";
    return;
  }
 
  if (dir != 'x' && dir != 'X' && dir != 'z' && dir != 'Z') {
    std::cerr << m_className << "::AddStripOnPlaneY:\n"
              << "    Invalid direction (" << dir << ").\n"
              << "    Only strips in x or z direction are possible.\n";
    return;
  }
 
  if (fabs(smax - smin) < Small) {
    std::cerr << m_className << "::AddStripOnPlaneY:\n"
              << "    Strip width must be greater than zero.\n";
    return;
  }
 
  Strip newStrip;
  newStrip.type = label;
  newStrip.ind = -1;
  newStrip.smin = std::min(smin, smax);
  newStrip.smax = std::max(smin, smax);
  newStrip.gap = gap > Small ? gap : -1.;
 
  int iplane = 2;
  if (m_ynplan[3]) {
    const double d2 = fabs(m_coplan[2] - y);
    const double d3 = fabs(m_coplan[3] - y);
    if (d3 < d2) iplane = 3;
  }
 
  if (dir == 'x' || dir == 'X') {
    m_planes[iplane].strips1.push_back(std::move(newStrip));
  } else {
    m_planes[iplane].strips2.push_back(std::move(newStrip));
  }
}

Referenced by main().

AddTube()

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

Add a tube.

Definition at line 901 of file ComponentAnalyticField.cc.

                                                             {
  // Check if the provided parameters make sense.
  if (radius <= 0.0) {
    std::cerr << m_className << "::AddTube: Unphysical tube dimension.\n";
    return;
  }
 
  if (nEdges < 3 && nEdges != 0) {
    std::cerr << m_className << "::AddTube: Unphysical number of tube edges ("
              << nEdges << ")\n";
    return;
  }
 
  // If there is already a tube defined, print a warning message.
  if (m_tube) {
    std::cout << m_className << "::AddTube:\n"
              << "    Warning: Existing tube settings will be overwritten.\n";
  }
 
  // Set the coordinate system.
  m_tube = true;
  m_polar = false;
 
  // Set the tube parameters.
  m_cotube = radius;
  m_cotube2 = radius * radius;
  m_vttube = voltage;
 
  m_ntube = nEdges;
 
  m_planes[4].type = label;
  m_planes[4].ind = -1;
 
  // Force recalculation of the capacitance and signal matrices.
  m_cellset = false;
  m_sigset = false;
}

Referenced by GarfieldPhysics::InitializePhysics().

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

                                                                  {
  // Check if the provided parameters make sense.
  if (diameter <= 0.) {
    std::cerr << m_className << "::AddWire: Unphysical wire diameter.\n";
    return;
  }
 
  if (tension <= 0.) {
    std::cerr << m_className << "::AddWire: Unphysical wire tension.\n";
    return;
  }
 
  if (rho <= 0.) {
    std::cerr << m_className << "::AddWire: Unphysical wire density.\n";
    return;
  }
 
  if (length <= 0.) {
    std::cerr << m_className << "::AddWire: Unphysical wire length.\n";
    return;
  }
 
  if (ntrap <= 0) {
    std::cerr << m_className << "::AddWire: Nbr. of trap radii must be > 0.\n";
    return;
  }
 
  if (m_polar && x <= diameter) {
    std::cerr << m_className << "::AddWire: Wire is too close to the origin.\n";
    return;
  }
 
  // Create a new wire.
  Wire wire;
  if (m_polar) {
    double r = 0., phi = 0.;
    Polar2Internal(x, y, r, phi);
    wire.x = r;
    wire.y = phi;
    wire.r = 0.5 * diameter / x; 
  } else {
    wire.x = x;
    wire.y = y;
    wire.r = 0.5 * diameter;
  }
  wire.v = voltage;
  wire.u = length;
  wire.type = label;
  wire.e = 0.;
  wire.ind = -1;
  wire.nTrap = ntrap;
  wire.tension = tension;
  wire.density = rho;
  // Add the wire to the list
  m_w.push_back(std::move(wire));
  ++m_nWires;
 
  // Force recalculation of the capacitance and signal matrices.
  m_cellset = false;
  m_sigset = false;
}

Referenced by GarfieldPhysics::InitializePhysics().

ClearCharges()

void Garfield::ComponentAnalyticField::ClearCharges

(

)

Remove all point charges.

Definition at line 1616 of file ComponentAnalyticField.cc.

                                          {
  m_ch3d.clear();
  m_nTermBessel = 10;
  m_nTermPoly = 100;
}

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

Definition at line 41 of file ComponentAnalyticField.hh.

                                           {
    m = nullptr;
    // Calculate the field.
    status = Field(x, y, z, ex, ey, ez, v, true);
    // If the field is ok, get the medium.
    if (status == 0) {
      m = m_geometry ? m_geometry->GetMedium(x, y, z) : m_medium;
      if (!m) {
        status = -6;
      } else if (!m->IsDriftable()) {
        status = -5;
      }
    }
  }

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,z	coordinates [cm].
ex,ey,ez	components of the electric field [V/cm].
m	pointer to the medium at this location.
status	status flag

Status flags:

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

Implements Garfield::Component.

Definition at line 24 of file ComponentAnalyticField.hh.

                                                                               {
    m = nullptr;
    // Calculate the field.
    double v = 0.;
    status = Field(x, y, z, ex, ey, ez, v, false);
    // If the field is ok, get the medium.
    if (status == 0) {
      m = m_geometry ? m_geometry->GetMedium(x, y, z) : m_medium;
      if (!m) {
        status = -6;
      } else if (!m->IsDriftable()) {
        status = -5;
      }
    }
  }

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

                                                                         {
  //-----------------------------------------------------------------------
  //   FFIELD - Subroutine calculating the electric field at a given wire
  //            position, as if the wire itself were not there but with
  //            the presence of its mirror images.
  //   VARIABLES : IW         : wire number
  //               EX, EY     : x- and y-component of the electric field.
  //   (Last changed on 27/ 1/96.)
  //-----------------------------------------------------------------------
  ex = ey = 0.;
  // Check the wire number.
  if (iw >= m_nWires) {
    std::cerr << m_className << "::ElectricFieldAtWire: Index out of range.\n";
    return false;
  }
  // Set the flags appropriately.
  std::vector<bool> cnalso(m_nWires, true);
  cnalso[iw] = false;
 
  const double xpos = m_w[iw].x;
  const double ypos = m_w[iw].y;
  // Call the appropriate function.
  switch (m_cellType) {
    case A00:
      FieldAtWireA00(xpos, ypos, ex, ey, cnalso);
      break;
    case B1X:
      FieldAtWireB1X(xpos, ypos, ex, ey, cnalso);
      break;
    case B1Y:
      FieldAtWireB1Y(xpos, ypos, ex, ey, cnalso);
      break;
    case B2X:
      FieldAtWireB2X(xpos, ypos, ex, ey, cnalso);
      break;
    case B2Y:
      FieldAtWireB2Y(xpos, ypos, ex, ey, cnalso);
      break;
    case C10:
      FieldAtWireC10(xpos, ypos, ex, ey, cnalso);
      break;
    case C2X:
      FieldAtWireC2X(xpos, ypos, ex, ey, cnalso);
      break;
    case C2Y:
      FieldAtWireC2Y(xpos, ypos, ex, ey, cnalso);
      break;
    case C30:
      FieldAtWireC30(xpos, ypos, ex, ey, cnalso);
      break;
    case D10:
      FieldAtWireD10(xpos, ypos, ex, ey, cnalso);
      break;
    case D20:
      FieldAtWireD20(xpos, ypos, ex, ey, cnalso);
      break;
    case D30:
      FieldAtWireD30(xpos, ypos, ex, ey, cnalso);
      break;
    default:
      std::cerr << m_className << "::ElectricFieldAtWire:\n"
                << "    Unknown cell type (id " << m_cellType << ")\n";
      return false;
  }
  // Correct for the equipotential planes.
  ex -= m_corvta;
  ey -= m_corvtb;
  if (m_polar) {
    const double r = exp(xpos);
    const double er = ex / r;
    const double ep = ey / r;
    const double ct = cos(ypos);
    const double st = sin(ypos);
    ex = +ct * er - st * ep;
    ey = +st * er + ct * ep;
  }
  return true;
}

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

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

                                                            {
 
  m_cellset = false;
  m_sigset = false;
  m_dipole = on;
}

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

{ 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

iw	index of the wire
xMap	coordinates of the grid lines in x
yMap	coordinates of the grid lines in y
fxMap	x-components of the force at the grid points
fyMap	y-components of the force at the grid points

Definition at line 1918 of file ComponentAnalyticField.cc.

                                          {
  if (!m_cellset && !Prepare()) {
    std::cerr << m_className << "::ForcesOnWire: Cell not set up.\n";
    return false;
  }
 
  if (iw >= m_nWires) {
    std::cerr << m_className << "::ForcesOnWire: Wire index out of range.\n";
    return false;
  }
  if (m_polar) {
    std::cerr << m_className << "::ForcesOnWire: Cannot handle polar cells.\n";
    return false;
  }
  const auto& wire = m_w[iw];
  // Compute a 'safe box' around the wire.
  double bxmin = m_perx ? wire.x - 0.5 * m_sx : 2 * m_xmin - m_xmax;
  double bxmax = m_perx ? wire.x + 0.5 * m_sx : 2 * m_xmax - m_xmin;
  double bymin = m_pery ? wire.y - 0.5 * m_sy : 2 * m_ymin - m_ymax;
  double bymax = m_pery ? wire.y + 0.5 * m_sy : 2 * m_ymax - m_ymin;
 
  // If the initial area is almost zero in 1 direction, make it square.
  if (std::abs(bxmax - bxmin) < 0.1 * std::abs(bymax - bymin)) {
    bxmin = wire.x - 0.5 * std::abs(bymax - bymin);
    bxmax = wire.x + 0.5 * std::abs(bymax - bymin);
  } else if (std::abs(bymax - bymin) < 0.1 * std::abs(bxmax - bxmin)) {
    bymin = wire.y - 0.5 * std::abs(bxmax - bxmin);
    bymax = wire.y + 0.5 * std::abs(bxmax - bxmin);
  }
  const double dw = 2 * wire.r;
  // Scan the other wires.
  for (unsigned int j = 0; j < m_nWires; ++j) {
    if (j == iw) continue;
    const double xj = m_w[j].x;
    const double yj = m_w[j].y;
    const double dj = 2 * m_w[j].r;
    double xnear = m_perx ? xj - m_sx * int(round((xj - wire.x) / m_sx)) : xj;
    double ynear = m_pery ? yj - m_sy * int(round((yj - wire.y) / m_sy)) : yj;
    if (std::abs(xnear - wire.x) > std::abs(ynear - wire.y)) {
      if (xnear < wire.x) {
        bxmin = std::max(bxmin, xnear + dj + dw);
        if (m_perx) bxmax = std::min(bxmax, xnear + m_sx - dj - dw);
      } else {
        bxmax = std::min(bxmax, xnear - dj - dw);
        if (m_perx) bxmin = std::max(bxmin, xnear - m_sx + dj + dw);
      }
    } else {
      if (ynear < wire.y) {
        bymin = std::max({bymin, ynear - dj - dw, ynear + dj + dw});
        if (m_pery) bymax = std::min(bymax, ynear + m_sy - dj - dw);
      } else {
        bymax = std::min({bymax, ynear - dj - dw, ynear + dj + dw});
        if (m_pery) bymin = std::max(bymin, ynear - m_sy + dj + dw);
      }
    }
  }
  // Scan the planes.
  if (m_ynplan[0]) bxmin = std::max(bxmin, m_coplan[0] + dw);
  if (m_ynplan[1]) bxmax = std::min(bxmax, m_coplan[1] - dw);
  if (m_ynplan[2]) bymin = std::max(bymin, m_coplan[2] + dw);
  if (m_ynplan[3]) bymax = std::min(bymax, m_coplan[3] - dw);
 
  // If there is a tube, check all corners.
  if (m_tube) {
    const double d2 = m_cotube2 - dw * dw;
    if (d2 < Small) {
      std::cerr << m_className << "::ForcesOnWire:\n    Diameter of wire " << iw
                << " is too large compared to the tube.\n";
      return false;
    }
 
    double corr = sqrt((bxmin * bxmin + bymin * bymin) / d2);
    if (corr > 1.) {
      bxmin /= corr;
      bymin /= corr;
    }
    corr = sqrt((bxmin * bxmin + bymax * bymax) / d2);
    if (corr > 1.) {
      bxmin /= corr;
      bymax /= corr;
    }
    corr = sqrt((bxmax * bxmax + bymin * bymin) / d2);
    if (corr > 1.) {
      bxmax /= corr;
      bymin /= corr;
    }
    corr = sqrt((bxmax * bxmax + bymax * bymax) / d2);
    if (corr > 1) {
      bxmax /= corr;
      bymax /= corr;
    }
  }
  // Make sure we found a reasonable 'safe area'.
  if ((bxmin - wire.x) * (wire.x - bxmax) <= 0 ||
      (bymin - wire.y) * (wire.y - bymax) <= 0) {
    std::cerr << m_className << "::ForcesOnWire:\n    Unable to find an area "
              << "free of elements around wire " << iw << ".\n";
    return false;
  }
  // Now set a reasonable scanning range.
  double sxmin = bxmin;
  double sxmax = bxmax;
  double symin = bymin;
  double symax = bymax;
  if (m_scanRange == ScanningRange::User) {
    // User-specified range.
    sxmin = wire.x + m_xScanMin;
    symin = wire.y + m_yScanMin;
    sxmax = wire.x + m_xScanMax;
    symax = wire.y + m_yScanMax;
  } else if (m_scanRange == ScanningRange::FirstOrder) {
    // Get the field and force at the nominal position.
    double ex = 0., ey = 0.;
    ElectricFieldAtWire(iw, ex, ey);
    double fx = -ex * wire.e * Internal2Newton;
    double fy = -ey * wire.e * Internal2Newton;
    if (m_useGravitationalForce) {
      // Mass per unit length [kg / cm].
      const double m = 1.e-3 * wire.density * Pi * wire.r * wire.r;
      fx -= m_down[0] * GravitationalAcceleration * m;
      fy -= m_down[1] * GravitationalAcceleration * m;
    }
    const double u2 = wire.u * wire.u;
    const double shiftx =
        -125 * fx * u2 / (GravitationalAcceleration * wire.tension);
    const double shifty =
        -125 * fy * u2 / (GravitationalAcceleration * wire.tension);
    // If 0th order estimate of shift is not small.
    const double tol = 0.1 * wire.r;
    if (std::abs(shiftx) > tol || std::abs(shifty) > tol) {
      sxmin = std::max(bxmin, std::min(wire.x + m_scaleRange * shiftx,
                                       wire.x - shiftx / m_scaleRange));
      symin = std::max(bymin, std::min(wire.y + m_scaleRange * shifty,
                                       wire.y - shifty / m_scaleRange));
      sxmax = std::min(bxmax, std::max(wire.x + m_scaleRange * shiftx,
                                       wire.x - shiftx / m_scaleRange));
      symax = std::min(bymax, std::max(wire.y + m_scaleRange * shifty,
                                       wire.y - shifty / m_scaleRange));
      // If one is very small, make the area square within bounds.
      if (std::abs(sxmax - sxmin) < 0.1 * std::abs(symax - symin)) {
        sxmin = std::max(bxmin, wire.x - 0.5 * std::abs(symax - symin));
        sxmax = std::min(bxmax, wire.x + 0.5 * std::abs(symax - symin));
      } else if (std::abs(symax - symin) < 0.1 * std::abs(sxmax - sxmin)) {
        symin = std::max(bymin, wire.y - 0.5 * std::abs(sxmax - sxmin));
        symax = std::min(bymax, wire.y + 0.5 * std::abs(sxmax - sxmin));
      }
    }
  }
  if (m_debug) {
    std::cout << m_className << "::ForcesOnWire:\n";
    std::printf("    Free area %12.5e < x < %12.5e\n", bxmin, bxmax);
    std::printf("              %12.5e < y < %12.5e\n", bymin, bymax);
    std::printf("    Scan area %12.5e < x < %12.5e\n", sxmin, sxmax);
    std::printf("              %12.5e < y < %12.5e\n", symin, symax);
  }
 
  xMap.resize(m_nScanX);
  const double stepx = (sxmax - sxmin) / (m_nScanX - 1);
  for (unsigned int i = 0; i < m_nScanX; ++i) {
    xMap[i] = sxmin + i * stepx;
  }
  yMap.resize(m_nScanY);
  const double stepy = (symax - symin) / (m_nScanY - 1);
  for (unsigned int i = 0; i < m_nScanY; ++i) {
    yMap[i] = symin + i * stepy;
  }
  // Save the original coordinates of the wire.
  const double x0 = wire.x;
  const double y0 = wire.y;
  // Prepare interpolation tables.
  fxMap.assign(m_nScanX, std::vector<double>(m_nScanY, 0.));
  fyMap.assign(m_nScanX, std::vector<double>(m_nScanY, 0.));
  bool ok = true;
  for (unsigned int i = 0; i < m_nScanX; ++i) {
    for (unsigned int j = 0; j < m_nScanY; ++j) {
      // Get the wire position for this shift.
      m_w[iw].x = xMap[i];
      m_w[iw].y = yMap[j];
      // Verify the current situation.
      if (!WireCheck()) {
        std::cerr << m_className << "::ForcesOnWire: Wire " << iw << ".\n"
                  << "    Scan involves a disallowed wire position.\n";
        ok = false;
        continue;
      }
      // Recompute the charges for this configuration.
      if (!Setup()) {
        std::cerr << m_className << "::ForcesOnWire: Wire " << iw << ".\n"
                  << "    Failed to compute charges at a scan point.\n";
        ok = false;
        continue;
      }
      // Compute the forces.
      double ex = 0., ey = 0.;
      ElectricFieldAtWire(iw, ex, ey);
      fxMap[i][j] = -ex * wire.e * Internal2Newton;
      fyMap[i][j] = -ey * wire.e * Internal2Newton;
    }
  }
  // Place the wire back in its shifted position.
  m_w[iw].x = x0;
  m_w[iw].y = y0;
  Setup();
  return ok;
}

Referenced by WireDisplacement().

GetBoundingBox()

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

overridevirtual

Get the bounding box coordinates.

Reimplemented from Garfield::Component.

Definition at line 319 of file ComponentAnalyticField.cc.

                                        {
  // If a geometry is present, try to get the bounding box from there.
  if (m_geometry) {
    if (m_geometry->GetBoundingBox(x0, y0, z0, x1, y1, z1)) return true;
  }
  // Otherwise, return the cell dimensions.
  return GetElementaryCell(x0, y0, z0, x1, y1, z1);
}

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

                          {
    if (!m_cellset) {
      if (CellCheck()) CellType();
    }
    return GetCellType(m_cellType);
  }

Referenced by GetCellType(), and PrintCell().

GetElementaryCell()

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

overridevirtual

Get the coordinates of the elementary cell.

Reimplemented from Garfield::Component.

Definition at line 330 of file ComponentAnalyticField.cc.

                                        {
  if (!m_cellset && !Prepare()) return false;
  if (m_polar) {
    double rmax, thetamax;
    Internal2Polar(m_xmax, m_ymax, rmax, thetamax);
    x0 = -rmax;
    y0 = -rmax;
    x1 = +rmax;
    y1 = +rmax;
  } else {
    x0 = m_xmin;
    y0 = m_ymin;
    x1 = m_xmax;
    y1 = m_ymax;
  }
  z0 = m_zmin;
  z1 = m_zmax;
  return true;
}

Referenced by GetBoundingBox().

GetGravity()

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

Get the gravity orientation.

Definition at line 1910 of file ComponentAnalyticField.cc.

                                                          {
 
  dx = m_down[0];
  dy = m_down[1];
  dz = m_down[2];
}

GetMedium()

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

overridevirtual

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

Reimplemented from Garfield::Component.

Definition at line 253 of file ComponentAnalyticField.cc.

                                                            {
  if (m_geometry) return m_geometry->GetMedium(xin, yin, zin);
 
  // Make sure the cell is prepared.
  if (!m_cellset && !Prepare()) return nullptr;
 
  double xpos = xin, ypos = yin;
  if (m_polar) Cartesian2Internal(xin, yin, xpos, ypos);
  // In case of periodicity, move the point into the basic cell.
  if (m_perx) {
    xpos -= m_sx * int(round(xin / m_sx));
  }
  double arot = 0.;
  if (m_pery && m_tube) {
    Cartesian2Polar(xin, yin, xpos, ypos);
    arot = RadToDegree * m_sy * int(round(DegreeToRad * ypos / m_sy));
    ypos -= arot;
    Polar2Cartesian(xpos, ypos, xpos, ypos);
  } else if (m_pery) {
    ypos -= m_sy * int(round(ypos / m_sy));
  }
 
  // Move the point to the correct side of the plane.
  if (m_perx && m_ynplan[0] && xpos <= m_coplan[0]) xpos += m_sx;
  if (m_perx && m_ynplan[1] && xpos >= m_coplan[1]) xpos -= m_sx;
  if (m_pery && m_ynplan[2] && ypos <= m_coplan[2]) ypos += m_sy;
  if (m_pery && m_ynplan[3] && ypos >= m_coplan[3]) ypos -= m_sy;
 
  // In case (xpos, ypos) is located behind a plane there is no field.
  if (m_tube) {
    if (!InTube(xpos, ypos, m_cotube, m_ntube)) return nullptr;
  } else {
    if ((m_ynplan[0] && xpos < m_coplan[0]) ||
        (m_ynplan[1] && xpos > m_coplan[1]) ||
        (m_ynplan[2] && ypos < m_coplan[2]) ||
        (m_ynplan[3] && ypos > m_coplan[3])) {
      return nullptr;
    }
  }
 
  // If (xpos, ypos) is within a wire, there is no field either.
  for (const auto& wire : m_w) {
    double dx = xpos - wire.x;
    double dy = ypos - wire.y;
    // Correct for periodicities.
    if (m_perx) dx -= m_sx * int(round(dx / m_sx));
    if (m_pery) dy -= m_sy * int(round(dy / m_sy));
    // Check the actual position.
    if (dx * dx + dy * dy < wire.r * wire.r) return nullptr;
  }
  return m_medium;
}

GetNumberOfPlanesPhi()

unsigned int Garfield::ComponentAnalyticField::GetNumberOfPlanesPhi

(

)

const

Get the number of equipotential planes at constant phi.

Definition at line 1669 of file ComponentAnalyticField.cc.

                                                                {
  if (!m_polar) {
    return 0;
  } else if (m_ynplan[2] && m_ynplan[3]) {
    return 2;
  } else if (m_ynplan[2] || m_ynplan[3]) {
    return 1;
  }
  return 0;
}

GetNumberOfPlanesR()

unsigned int Garfield::ComponentAnalyticField::GetNumberOfPlanesR

(

)

const

Get the number of equipotential planes at constant radius.

Definition at line 1658 of file ComponentAnalyticField.cc.

                                                              {
  if (!m_polar) {
    return 0;
  } else if (m_ynplan[0] && m_ynplan[1]) {
    return 2;
  } else if (m_ynplan[0] || m_ynplan[1]) {
    return 1;
  }
  return 0;
}

GetNumberOfPlanesX()

unsigned int Garfield::ComponentAnalyticField::GetNumberOfPlanesX

(

)

const

Get the number of equipotential planes at constant x.

Definition at line 1636 of file ComponentAnalyticField.cc.

                                                              {
  if (m_polar) {
    return 0;
  } else if (m_ynplan[0] && m_ynplan[1]) {
    return 2;
  } else if (m_ynplan[0] || m_ynplan[1]) {
    return 1;
  }
  return 0;
}

GetNumberOfPlanesY()

unsigned int Garfield::ComponentAnalyticField::GetNumberOfPlanesY

(

)

const

Get the number of equipotential planes at constant y.

Definition at line 1647 of file ComponentAnalyticField.cc.

                                                              {
  if (m_polar) {
    return 0;
  } else if (m_ynplan[2] && m_ynplan[3]) {
    return 2;
  } else if (m_ynplan[2] || m_ynplan[3]) {
    return 1;
  }
  return 0;
}

GetNumberOfWires()

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

inline

Get the number of wires.

Definition at line 243 of file ComponentAnalyticField.hh.

{ return m_nWires; }

GetPeriodicityPhi()

bool Garfield::ComponentAnalyticField::GetPeriodicityPhi

(

double &

s

)

Get the periodicity [degree] in phi.

Definition at line 1517 of file ComponentAnalyticField.cc.

                                                        {
  if (!m_periodic[1] || (!m_polar && !m_tube)) {
    s = 0.;
    return false;
  }
 
  s = RadToDegree * m_sy;
  return true;
}

GetPeriodicityX()

bool Garfield::ComponentAnalyticField::GetPeriodicityX

(

double &

s

)

Get the periodic length in the x-direction.

Definition at line 1497 of file ComponentAnalyticField.cc.

                                                      {
  if (!m_periodic[0] || m_polar) {
    s = 0.;
    return false;
  }
 
  s = m_sx;
  return true;
}

GetPeriodicityY()

bool Garfield::ComponentAnalyticField::GetPeriodicityY

(

double &

s

)

Get the periodic length in the y-direction.

Definition at line 1507 of file ComponentAnalyticField.cc.

                                                      {
  if (!m_periodic[1] || m_polar) {
    s = 0.;
    return false;
  }
 
  s = m_sy;
  return true;
}

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

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

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

                                                               {
  if (!m_polar || i >= 2 || (i == 1 && !m_ynplan[1])) {
    std::cerr << m_className << "::GetPlaneR: Index out of range.\n";
    return false;
  }
 
  r = exp(m_coplan[i]);
  voltage = m_vtplan[i];
  label = m_planes[i].type;
  return true;
}

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

                                                               {
  if (m_polar || i >= 2 || (i == 1 && !m_ynplan[1])) {
    std::cerr << m_className << "::GetPlaneX: Index out of range.\n";
    return false;
  }
 
  x = m_coplan[i];
  voltage = m_vtplan[i];
  label = m_planes[i].type;
  return true;
}

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

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

GetTube()

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

Retrieve the tube parameters.

Definition at line 1764 of file ComponentAnalyticField.cc.

                                                             {
  if (!m_tube) return false;
  r = m_cotube;
  voltage = m_vttube;
  nEdges = m_ntube;
  label = m_planes[4].type;
  return true;
}

GetVoltageRange()

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

overridevirtual

Calculate the voltage range [V].

Implements Garfield::Component.

Definition at line 307 of file ComponentAnalyticField.cc.

                                                                       {
  // Make sure the cell is prepared.
  if (!m_cellset && !Prepare()) {
    std::cerr << m_className << "::GetVoltageRange: Cell not set up.\n";
    return false;
  }
 
  pmin = m_vmin;
  pmax = m_vmax;
  return true;
}

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

                                                                       {
  if (i >= m_nWires) {
    std::cerr << m_className << "::GetWire: Index out of range.\n";
    return false;
  }
 
  if (m_polar) {
    double r = 0., theta = 0.;
    Internal2Polar(m_w[i].x, m_w[i].y, r, theta);
    x = r;
    y = theta;
    diameter = 2 * m_w[i].r * r;
  } else {
    x = m_w[i].x;
    y = m_w[i].y;
    diameter = 2 * m_w[i].r;
  }
  voltage = m_w[i].v;
  label = m_w[i].type;
  length = m_w[i].u;
  charge = m_w[i].e;
  ntrap = m_w[i].nTrap;
  return true;
}

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

q0	charge of the particle [in elementary charges].
x0,y0,z0	position [cm] of the particle.
xw,yw	coordinates of the wire (if applicable).
rw	radius of the wire (if applicable).

Reimplemented from Garfield::Component.

Definition at line 751 of file ComponentAnalyticField.cc.

                                                        {
 
  double x0 = xin;
  double y0 = yin;
  if (m_polar) {
    Cartesian2Internal(xin, yin, x0, y0);
  }
  // In case of periodicity, move the point into the basic cell.
  int nX = 0, nY = 0, nPhi = 0;
  if (m_perx) {
    nX = int(round(x0 / m_sx));
    x0 -= m_sx * nX;
  }
  if (m_pery && m_tube) {
    Cartesian2Polar(xin, yin, x0, y0);
    nPhi = int(round(DegreeToRad * y0 / m_sy));
    y0 -= RadToDegree * m_sy * nPhi;
    Polar2Cartesian(x0, y0, x0, y0);
  } else if (m_pery) {
    nY = int(round(y0 / m_sy));
    y0 -= m_sy * nY;
  }
  // Move the point to the correct side of the plane.
  std::array<bool, 4> shift = {false, false, false, false};
  if (m_perx && m_ynplan[0] && x0 <= m_coplan[0]) {
    x0 += m_sx;
    shift[0] = true;
  }
  if (m_perx && m_ynplan[1] && x0 >= m_coplan[1]) {
    x0 -= m_sx;
    shift[1] = true;
  }
  if (m_pery && m_ynplan[2] && y0 <= m_coplan[2]) {
    y0 += m_sy;
    shift[2] = true;
  }
  if (m_pery && m_ynplan[3] && y0 >= m_coplan[3]) {
    y0 -= m_sy;
    shift[3] = true;
  }
 
  for (const auto& wire : m_w) {
    // Skip wires with the wrong charge.
    if (qin * wire.e > 0.) continue;
    const double dxw0 = wire.x - x0;
    const double dyw0 = wire.y - y0;
    const double r2 = dxw0 * dxw0 + dyw0 * dyw0;
    const double rTrap = wire.r * wire.nTrap;
    if (r2 < rTrap * rTrap) {
      xw = wire.x;
      yw = wire.y;
      rw = wire.r;
      if (shift[0]) xw -= m_sx;
      if (shift[1]) xw += m_sx;
      if (shift[2]) yw -= m_sy;
      if (shift[3]) yw += m_sy;
      if (m_pery && m_tube) {
        double rhow, phiw;
        Cartesian2Polar(xw, yw, rhow, phiw);
        phiw += RadToDegree * m_sy * nPhi;
        Polar2Cartesian(rhow, phiw, xw, yw);
      } else if (m_pery) {
        y0 += m_sy * nY;
      }
      if (m_perx) xw += m_sx * nX;
      if (m_polar) {
        Internal2Cartesian(xw, yw, xw, yw);
        rw *= exp(wire.x);
      }
      if (m_debug) {
        std::cout << m_className << "::IsInTrapRadius: (" << xin << ", "
                  << yin << ", " << zin << ")" << " within trap radius.\n";
      }
      return true;
    }
  }
 
  return false;
}

IsPolar()

bool Garfield::ComponentAnalyticField::IsPolar ( ) const

inline

Are polar coordinates being used?

Definition at line 181 of file ComponentAnalyticField.hh.

{ 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,z0	first point [cm].
x1,y1,z1	second point [cm]
xc,yc,zc	point [cm] where the line crosses the wire or the coordinates of the wire centre.
centre	flag whether to return the coordinates of the line-wire crossing point or of the wire centre.
rc	radius [cm] of the wire.

Reimplemented from Garfield::Component.

Definition at line 650 of file ComponentAnalyticField.cc.

                                                                          {
  xc = xx0;
  yc = yy0;
  zc = z0;
 
  if (m_w.empty()) return false;
 
  double x0 = xx0;
  double y0 = yy0;
  double x1 = xx1;
  double y1 = yy1;
  if (m_polar) {
    Cartesian2Internal(xx0, yy0, x0, y0);
    Cartesian2Internal(xx1, yy1, x1, y1);
  }
  const double dx = x1 - x0;
  const double dy = y1 - y0;
  const double d2 = dx * dx + dy * dy;
  // Check that the step length is non-zero.
  if (d2 < Small) return false;
  const double invd2 = 1. / d2;
 
  // Check if a whole period has been crossed.
  if ((m_perx && fabs(dx) >= m_sx) || (m_pery && fabs(dy) >= m_sy)) {
    std::cerr << m_className << "::IsWireCrossed:\n"
              << "    Particle crossed more than one period.\n";
    return false;
  }
 
  // Both coordinates are assumed to be located inside
  // the drift area and inside a drift medium.
  // This should have been checked before this call.
 
  const double xm = 0.5 * (x0 + x1);
  const double ym = 0.5 * (y0 + y1);
  double dMin2 = 0.;
  for (const auto& wire : m_w) {
    double xw = wire.x;
    if (m_perx) {
      xw += m_sx * int(round((xm - xw) / m_sx));
    }
    double yw = wire.y;
    if (m_pery) {
      yw += m_sy * int(round((ym - yw) / m_sy));
    }
    // Calculate the smallest distance between track and wire.
    const double xIn0 = dx * (xw - x0) + dy * (yw - y0);
    // Check if the minimum is located before (x0, y0).
    if (xIn0 < 0.) continue;
    const double xIn1 = -(dx * (xw - x1) + dy * (yw - y1));
    // Check if the minimum is located behind (x1, y1).
    if (xIn1 < 0.) continue;
    // Minimum is located between (x0, y0) and (x1, y1).
    const double xw0 = xw - x0;
    const double xw1 = xw - x1;
    const double yw0 = yw - y0;
    const double yw1 = yw - y1;
    const double dw02 = xw0 * xw0 + yw0 * yw0;
    const double dw12 = xw1 * xw1 + yw1 * yw1;
    if (xIn1 * xIn1 * dw02 > xIn0 * xIn0 * dw12) {
      dMin2 = dw02 - xIn0 * xIn0 * invd2;
    } else {
      dMin2 = dw12 - xIn1 * xIn1 * invd2;
    }
    // Add in the times nTrap to account for the trap radius.
    const double r2 = wire.r * wire.r;
    if (dMin2 < r2) {
      // Wire has been crossed.
      if (centre) {
        if (m_polar) {
          Internal2Cartesian(xw, yw, xc, yc);
        } else {
          xc = xw;
          yc = yw;
        }
      } else {
        // Find the point of intersection.
        const double p = -xIn0 * invd2;
        const double q = (dw02 - r2) * invd2;
        const double s = sqrt(p * p - q);
        const double t = std::min(-p + s, -p - s);
        if (m_polar) {
          Internal2Cartesian(x0 + t * dx, y0 + t * dy, xc, yc);
        } else {
          xc = x0 + t * dx;
          yc = y0 + t * dy;
        }
        zc = z0 + t * (z1 - z0);
      }
      rc = wire.r;
      if (m_polar) rc *= exp(wire.x);
      return true;
    }
  }
  return false;
}

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

iw	Index of the wire.
order	Order of the highest multipole moment.
print	Print information about the fitting process.
plot	Plot the potential and multipole fit around the wire.
rmult	Distance in multiples of the wire radius at which the decomposition is to be carried out.
eps	Used in the fit for calculating the covariance matrix.
nMaxIter	Maximum number of iterations in the fit.

Definition at line 10400 of file ComponentAnalyticField.cc.

                                 {
 
  //-----------------------------------------------------------------------
  //   EFMWIR - Computes the dipole moment of a given wire.
  //-----------------------------------------------------------------------
  if (!m_cellset && !Prepare()) return false;
  // Check input parameters.
  if (iw >= m_nWires) {
    std::cerr << m_className << "::MultipoleMoments:\n"
              << "    Wire index out of range.\n";
    return false;
  }
  if (eps <= 0.) {
    std::cerr << m_className << "::MultipoleMoments:\n"
              << "    Epsilon must be positive.\n";
    return false;
  }
  if (nPoles == 0) {
    std::cerr << m_className << "::MultipoleMoments:\n"
              << "    Multipole order out of range.\n";
    return false;
  }
  if (rmult <= 0.) {
    std::cerr << m_className << "::MultipoleMoments:\n"
              << "    Radius multiplication factor out of range.\n";
    return false;
  }
 
  const double xw = m_w[iw].x;
  const double yw = m_w[iw].y;
  // Set the radius of the wire to 0.
  const double rw = m_w[iw].r;
  m_w[iw].r = 0.;
 
  // Loop around the wire.
  constexpr unsigned int nPoints = 20000;
  std::vector<double> angle(nPoints, 0.);
  std::vector<double> volt(nPoints, 0.);
  std::vector<double> weight(nPoints, 1.);
  for (unsigned int i = 0; i < nPoints; ++i) {
    // Set angle around wire.
    angle[i] = TwoPi * (i + 1.) / nPoints;
    // Compute E field, make sure the point is in a free region.
    const double x = xw + rmult * rw * cos(angle[i]); 
    const double y = yw + rmult * rw * sin(angle[i]);
    double ex = 0., ey = 0., ez = 0., v = 0.;
    if (Field(x, y, 0., ex, ey, ez, v, true) != 0) { 
      std::cerr << m_className << "::MultipoleMoments:\n"
                << "    Unexpected location code. Computation stopped.\n";
      m_w[iw].r = rw;
      return false;
    }
    // Assign the result to the fitting array.
    volt[i] = v;
  } 
  // Restore the wire diameter.
  m_w[iw].r = rw;
 
  // Determine the maximum, minimum and average.
  double vmin = *std::min_element(volt.cbegin(), volt.cend());
  double vmax = *std::max_element(volt.cbegin(), volt.cend());
  double vave = std::accumulate(volt.cbegin(), volt.cend(), 0.) / nPoints;
  // Subtract the wire potential to centre the data more or less.
  for (unsigned int i = 0; i < nPoints; ++i) volt[i] -= vave;
  vmax -= vave;
  vmin -= vave;
 
  // Perform the fit.
  const double vm = 0.5 * fabs(vmin) + fabs(vmax);
  double chi2 = 1.e-6 * nPoints * vm * vm;
  const double dist = 1.e-3 * (1. + vm);
  const unsigned int nPar = 2 * nPoles + 1;
  std::vector<double> pars(nPar, 0.);
  std::vector<double> epar(nPar, 0.);
  pars[0] = 0.5 * (vmax + vmin);
  for (unsigned int i = 1; i <= nPoles; ++i) {
    pars[2 * i - 1] = 0.5 * (vmax - vmin);
    pars[2 * i] = 0.;
  } 
 
  auto f = [nPoles](const double x, const std::vector<double>& par) {
    // EFMFUN
    // Sum the series, initial value is the monopole term.
    double sum = par[0];
    for (unsigned int k = 1; k <= nPoles; ++k) {
      // Obtain the Legendre polynomial of this order and add to the series.
      const float cphi = cos(x - par[2 * k]);
      sum += par[2 * k - 1] * sqrt(k + 0.5) * Numerics::Legendre(k, cphi);
    }
    return sum;
  };
 
  if (!Numerics::LeastSquaresFit(f, pars, epar, angle, volt, weight, 
                                 nMaxIter, dist, chi2, eps, m_debug, print)) {
    std::cerr << m_className << "::MultipoleMoments:\n"
              << "    Fitting the multipoles failed; computation stopped.\n";
  }
  // Plot the result of the fit.
  if (plot) {
    const std::string name = ViewBase::FindUnusedCanvasName("cMultipole");
    TCanvas* cfit = new TCanvas(name.c_str(), "Multipoles");
    cfit->SetGridx();
    cfit->SetGridy();
    cfit->DrawFrame(0., vmin, TwoPi, vmax,
                    ";Angle around the wire [rad]; Potential - average [V]");
    TGraph graph;
    graph.SetLineWidth(2);
    graph.SetLineColor(kBlack);
    graph.DrawGraph(angle.size(), angle.data(), volt.data(), "lsame");
    // Sum of contributions.
    constexpr unsigned int nP = 1000;
    std::array<double, nP> xp;
    std::array<double, nP> yp;
    for (unsigned int i = 0; i < nP; ++i) {
      xp[i] = TwoPi * (i + 1.) / nP;
      yp[i] = f(xp[i], pars);
    }
    graph.SetLineColor(kViolet + 3);
    graph.DrawGraph(nP, xp.data(), yp.data(), "lsame");
    // Individual contributions.
    std::vector<double> parres = pars;
    for (unsigned int i = 1; i <= nPoles; ++i) parres[2 * i - 1] = 0.;
    for (unsigned int j = 1; j <= nPoles; ++j) {
      parres[2 * j - 1] = pars[2 * j - 1];
      for (unsigned int i = 0; i < nP; ++i) {
        yp[i] = f(xp[i], parres);
      }
      parres[2 * j - 1] = 0.;
      graph.SetLineColor(kAzure + j);
      graph.DrawGraph(nP, xp.data(), yp.data(), "lsame"); 
    }
    gPad->Update();
  }
 
  // Print the results.
  std::cout << m_className << "::MultipoleMoments:\n"
            << "    Multipole moments for wire " << iw << ":\n"
            << "  Moment            Value       Angle [degree]\n";
  std::printf("  %6u  %15.8f        Arbitrary\n", 0, vave);
  for (unsigned int i = 1; i <= nPoles; ++i) {
    // Remove radial term from the multipole moment.
    const double val = pow(rmult * rw, i) * pars[2 * i - 1];
    const double phi = RadToDegree * fmod(pars[2 * i], Pi);
    std::printf("  %6u  %15.8f  %15.8f\n", i, val, phi);
  }
  return true;
}

PrintCell()

void Garfield::ComponentAnalyticField::PrintCell

(

)

Print all available information on the cell.

Definition at line 352 of file ComponentAnalyticField.cc.

                                       {
  //-----------------------------------------------------------------------
  //   CELPRT - Subroutine printing all available information on the cell.
  //-----------------------------------------------------------------------
 
  // Make sure the cell is prepared.
  if (!m_cellset && !Prepare()) {
    std::cerr << m_className << "::PrintCell: Cell not set up.\n";
    return;
  }
  std::cout << m_className
            << "::PrintCell: Cell identification: " << GetCellType() << "\n";
  // Print positions of wires, applied voltages and resulting charges.
  if (!m_w.empty()) {
    std::cout << "  Table of the wires\n";
    if (m_polar) {
      std::cout << "  Nr    Diameter     r       phi     Voltage";
    } else {
      std::cout << "  Nr    Diameter     x        y      Voltage";
    }
    std::cout << "      Charge    Tension   Length   Density  Label\n";
    if (m_polar) {
      std::cout << "       [micron]     [cm]     [deg]    [Volt]";
    } else {
      std::cout << "       [micron]     [cm]     [cm]     [Volt]";
    }
    std::cout << "     [pC/cm]       [g]      [cm]    [g/cm3]\n";
    for (unsigned int i = 0; i < m_nWires; ++i) {
      const auto& w = m_w[i];
      double xw = w.x;
      double yw = w.y;
      double dw = 2 * w.r;
      if (m_polar) {
        Internal2Polar(w.x, w.y, xw, yw);
        dw *= xw;
      }
      std::printf(
          "%4d %9.2f %9.4f %9.4f %9.3f %12.4f %9.2f %9.2f %9.2f \"%s\"\n", i,
          1.e4 * dw, xw, yw, w.v, w.e * TwoPiEpsilon0 * 1.e-3, w.tension,
          w.u, w.density, w.type.c_str());
    }
  }
  // Print information on the tube if present.
  if (m_tube) {
    std::string shape;
    if (m_ntube == 0) {
      shape = "Circular";
    } else if (m_ntube == 3) {
      shape = "Triangular";
    } else if (m_ntube == 4) {
      shape = "Square";
    } else if (m_ntube == 5) {
      shape = "Pentagonal";
    } else if (m_ntube == 6) {
      shape = "Hexagonal";
    } else if (m_ntube == 7) {
      shape = "Heptagonal";
    } else if (m_ntube == 8) {
      shape = "Octagonal";
    } else {
      shape = "Polygonal with " + std::to_string(m_ntube) + " corners";
    }
    std::cout << "  Enclosing tube\n"
              << "    Potential:  " << m_vttube << " V\n"
              << "    Radius:     " << m_cotube << " cm\n"
              << "    Shape:      " << shape << "\n"
              << "    Label:      " << m_planes[4].type << "\n";
  }
  // Print data on the equipotential planes.
  if (m_ynplan[0] || m_ynplan[1] || m_ynplan[2] || m_ynplan[3]) {
    std::cout << "  Equipotential planes\n";
    // First those at const x or r.
    const std::string xr = m_polar ? "r" : "x";
    if (m_ynplan[0] && m_ynplan[1]) {
      std::cout << "    There are two planes at constant " << xr << ":\n";
    } else if (m_ynplan[0] || m_ynplan[1]) {
      std::cout << "    There is one plane at constant " << xr << ":\n";
    }
    for (unsigned int i = 0; i < 2; ++i) {
      if (!m_ynplan[i]) continue;
      if (m_polar) {
        std::cout << "    r = " << exp(m_coplan[i]) << " cm, ";
      } else {
        std::cout << "    x = " << m_coplan[i] << " cm, ";
      }
      if (fabs(m_vtplan[i]) > 1.e-4) {
        std::cout << "potential = " << m_vtplan[i] << " V, ";
      } else {
        std::cout << "earthed, ";
      }
      const auto& plane = m_planes[i];
      if (plane.type.empty() && plane.type != "?") {
        std::cout << "label = " << plane.type << ", ";
      }
      const unsigned int nStrips = plane.strips1.size() + plane.strips2.size();
      const unsigned int nPixels = plane.pixels.size();
      if (nStrips == 0 && nPixels == 0) {
        std::cout << "no strips or pixels.\n";
      } else if (nPixels == 0) {
        std::cout << nStrips << " strips.\n";
      } else if (nStrips == 0) {
        std::cout << nPixels << " pixels.\n";
      } else {
        std::cout << nStrips << " strips, " << nPixels << " pixels.\n";
      }
      for (const auto& strip : plane.strips2) {
        std::cout << "      ";
        if (m_polar) {
          double gap = i == 0 ? expm1(strip.gap) : -expm1(-strip.gap);
          gap *= exp(m_coplan[i]);
          std::cout << RadToDegree * strip.smin << " < phi < "
                    << RadToDegree * strip.smax
                    << " degrees, gap = " << gap << " cm";
        } else {
          std::cout << strip.smin << " < y < " << strip.smax
                    << " cm, gap = " << strip.gap << " cm";
        }
        if (!strip.type.empty() && strip.type != "?") {
          std::cout << " (\"" << strip.type << "\")";
        }
        std::cout << "\n";
      }
      for (const auto& strip : plane.strips1) {
        std::cout << "      " << strip.smin << " < z < " << strip.smax;
        if (m_polar) {
          double gap = i == 0 ? expm1(strip.gap) : -expm1(-strip.gap);
          gap *= exp(m_coplan[i]);
          std::cout << " cm, gap = " << gap << " cm";
        } else {
          std::cout << " cm, gap = " << strip.gap << " cm";
        }
        if (!strip.type.empty() && strip.type != "?") {
          std::cout << " (\"" << strip.type << "\")";
        }
        std::cout << "\n";
      }
      for (const auto& pix : plane.pixels) {
        std::cout << "      ";
        if (m_polar) {
          std::cout << RadToDegree * pix.smin << " < phi < "
                    << RadToDegree * pix.smax << " degrees, "; 
        } else {
          std::cout << pix.smin << " < y < " << pix.smax << " cm, ";
        }
        std::cout << pix.zmin << " < z < " << pix.zmax << " cm, gap = ";
        if (m_polar) {
          double gap = i == 0 ? expm1(pix.gap) : -expm1(-pix.gap);
          gap *= exp(m_coplan[i]);
          std::cout << gap << " cm";
        } else {
          std::cout << pix.gap << " cm";
        }
        if (!pix.type.empty() && pix.type != "?") {
          std::cout << " (\"" << pix.type << "\")";
        }
        std::cout << "\n";
      }
    }
    // Next the planes at constant y or phi
    const std::string yphi = m_polar ? "phi" : "y";
    if (m_ynplan[2] && m_ynplan[3]) {
      std::cout << "    There are two planes at constant " << yphi << ":\n";
    } else if (m_ynplan[2] || m_ynplan[3]) {
      std::cout << "    There is one plane at constant " << yphi << ":\n";
    }
    for (unsigned int i = 2; i < 4; ++i) {
      if (!m_ynplan[i]) continue;
      if (m_polar) {
        std::cout << "    phi = " << RadToDegree * m_coplan[i] << " degrees, ";
      } else {
        std::cout << "    y = " << m_coplan[i] << " cm, ";
      }
      if (fabs(m_vtplan[i]) > 1.e-4) {
        std::cout << "potential = " << m_vtplan[i] << " V, ";
      } else {
        std::cout << "earthed, ";
      }
      const auto& plane = m_planes[i];
      if (plane.type.empty() && plane.type != "?") {
        std::cout << "label = " << plane.type << ", ";
      }
      const unsigned int nStrips = plane.strips1.size() + plane.strips2.size();
      const unsigned int nPixels = plane.pixels.size();
      if (nStrips == 0 && nPixels == 0) {
        std::cout << "no strips or pixels.\n";
      } else if (nPixels == 0) {
        std::cout << nStrips << " strips.\n";
      } else if (nStrips == 0) {
        std::cout << nPixels << " pixels.\n";
      } else {
        std::cout << nStrips << " strips, " << nPixels << " pixels.\n";
      }
      for (const auto& strip : plane.strips2) {
        std::cout << "      ";
        if (m_polar) {
          std::cout << exp(strip.smin) << " < r < " << exp(strip.smax) 
                    << " cm, gap = " << RadToDegree * strip.gap << " degrees";
        } else {
          std::cout << strip.smin << " < x < " << strip.smax
                    << " cm, gap = " << strip.gap << " cm";
        }
        if (!strip.type.empty() && strip.type != "?") {
          std::cout << " (\"" << strip.type << "\")";
        }
        std::cout << "\n";
      }
      for (const auto& strip : plane.strips1) {
        std::cout << "      " << strip.smin << " < z < " << strip.smax;
        if (m_polar) {
          std::cout << " cm, gap = " << RadToDegree * strip.gap << " degrees";
        } else {
          std::cout << " cm, gap = " << strip.gap << " cm";
        }
        if (!strip.type.empty() && strip.type != "?") {
          std::cout << " (\"" << strip.type << "\")";
        }
        std::cout << "\n";
      }
      for (const auto& pix : plane.pixels) {
        std::cout << "      ";
        if (m_polar) {
          std::cout << exp(pix.smin) << " < r < " << exp(pix.smax) << " cm, "; 
        } else {
          std::cout << pix.smin << " < x < " << pix.smax << " cm, ";
        }
        std::cout << pix.zmin << " < z < " << pix.zmax << " cm, gap = ";
        if (m_polar) {
          std::cout << RadToDegree * pix.gap << " degrees";
        } else {
          std::cout << pix.gap << " cm";
        }
        if (!pix.type.empty() && pix.type != "?") {
          std::cout << " (\"" << pix.type << "\")";
        }
        std::cout << "\n";
      }
    }
  }
  // Print the type of periodicity.
  std::cout << "  Periodicity\n";
  if (m_perx) {
    std::cout << "    The cell is repeated every ";
    if (m_polar) {
      std::cout << exp(m_sx) << " cm in r.\n";
    } else {
      std::cout << m_sx << " cm in x.\n";
    }
  } else {
    if (m_polar) {
      std::cout << "    The cell is not periodic in r.\n";
    } else {
      std::cout << "    The cell has no translation periodicity in x.\n";
    }
  }
  if (m_pery) {
    std::cout << "    The cell is repeated every ";
    if (m_polar) {
      std::cout << RadToDegree * m_sy << " degrees in phi.\n";
    } else {
      std::cout << m_sy << " cm in y.\n";
    }
  } else {
    if (m_polar) {
      std::cout << "    The cell is not periodic in phi.\n";
    } else {
      std::cout << "    The cell has no translation periodicity in y.\n";
    }
  }
  std::cout << "  Other data\n";
  std::cout << "    Gravity vector: (" << m_down[0] << ", " << m_down[1]
            << ", " << m_down[2] << ") [g].\n";
  std::cout << "  Cell dimensions:\n    ";
  if (!m_polar) {
    std::cout << m_xmin << " < x < " << m_xmax << " cm, " << m_ymin << " < y < "
              << m_ymax << " cm.\n";
  } else {
    double xminp, yminp;
    Internal2Polar(m_xmin, m_ymin, xminp, yminp);
    double xmaxp, ymaxp;
    Internal2Polar(m_xmax, m_ymax, xmaxp, ymaxp);
    std::cout << xminp << " < r < " << xmaxp << " cm, " << yminp << " < phi < "
              << ymaxp << " degrees.\n";
  }
  std::cout << "  Potential range:\n    " << m_vmin << " < V < " << m_vmax
            << " V.\n";
  // Print voltage shift in case no equipotential planes are present.
  if (!(m_ynplan[0] || m_ynplan[1] || m_ynplan[2] || m_ynplan[3] || m_tube)) {
    std::cout << "  All voltages have been shifted by " << m_v0
              << " V to avoid net wire charge.\n";
  } else {
    // Print the net charge on wires.
    double sum = 0.;
    for (const auto& w : m_w) sum += w.e;
    std::cout << "  The net charge on the wires is "
              << 1.e-3 * TwoPiEpsilon0 * sum << " pC/cm.\n";
  }
}

PrintCharges()

void Garfield::ComponentAnalyticField::PrintCharges

(

)

const

Print a list of the point charges.

Definition at line 1622 of file ComponentAnalyticField.cc.

                                                {
  std::cout << m_className << "::PrintCharges:\n";
  if (m_ch3d.empty()) {
    std::cout << "    No charges present.\n";
    return;
  }
  std::cout << "      x [cm]      y [cm]      z [cm]      charge [fC]\n";
  for (const auto& charge : m_ch3d) {
    std::cout << "     " << std::setw(9) << charge.x << "   " << std::setw(9)
              << charge.y << "   " << std::setw(9) << charge.z << "   "
              << std::setw(11) << charge.e * FourPiEpsilon0 << "\n";
  }
}

SetCartesianCoordinates()

void Garfield::ComponentAnalyticField::SetCartesianCoordinates

(

)

Use Cartesian coordinates (default).

Definition at line 1584 of file ComponentAnalyticField.cc.

                                                     {
  if (m_polar) {
    std::cout << m_className << "::SetCartesianCoordinates:\n    "
              << "Switching to Cartesian coordinates; resetting the cell.\n";
    CellInit();
  }
  m_polar = false;
}

SetGravity()

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

Set the gravity orientation.

Definition at line 1896 of file ComponentAnalyticField.cc.

                                                         {
 
  const double d = sqrt(dx * dx + dy * dy + dz * dz);
  if (d > 0.) {
    m_down[0] = dx / d;
    m_down[1] = dy / d;
    m_down[2] = dz / d;
  } else {
    std::cerr << m_className << "::SetGravity:\n"
              << "    The gravity vector has zero norm ; ignored.\n";
  } 
} 

SetMedium()

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

inline

Set the medium inside the cell.

Definition at line 89 of file ComponentAnalyticField.hh.

{ m_medium = medium; }

Referenced by GarfieldPhysics::InitializePhysics(), and main().

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

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

                                                                  {
  if (n == 0) {
    std::cerr << m_className << "::SetNumberOfSteps:\n"
              << "    Number of steps must be > 0.\n";
    return;
  }
  m_nSteps = n;
}

SetPeriodicityPhi()

void Garfield::ComponentAnalyticField::SetPeriodicityPhi

(

const double

phi

)

Set the periodicity [degree] in phi.

Definition at line 1479 of file ComponentAnalyticField.cc.

                                                             {
  if (!m_polar && !m_tube) {
    std::cerr << m_className << "::SetPeriodicityPhi:\n"
              << "    Cannot use phi-periodicity with Cartesian coordinates.\n";
    return;
  }
  if (std::abs(360. - s * int(360. / s)) > 1.e-4) {
    std::cerr << m_className << "::SetPeriodicityPhi:\n"
              << "    Phi periods must divide 360; ignored.\n";
    return;
  }
 
  m_periodic[1] = true;
  m_sy = DegreeToRad * s;
  m_mtube = int(360. / s);
  UpdatePeriodicity();
}

SetPeriodicityX()

void Garfield::ComponentAnalyticField::SetPeriodicityX

(

const double

s

)

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

Definition at line 1445 of file ComponentAnalyticField.cc.

                                                           {
  if (m_polar) {
    std::cerr << m_className << "::SetPeriodicityX:\n"
              << "    Cannot use x-periodicity with polar coordinates.\n";
    return;
  }
  if (s < Small) {
    std::cerr << m_className << "::SetPeriodicityX:\n"
              << "    Periodic length must be greater than zero.\n";
    return;
  }
 
  m_periodic[0] = true;
  m_sx = s;
  UpdatePeriodicity();
}

SetPeriodicityY()

void Garfield::ComponentAnalyticField::SetPeriodicityY

(

const double

s

)

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

Definition at line 1462 of file ComponentAnalyticField.cc.

                                                           {
  if (m_polar) {
    std::cerr << m_className << "::SetPeriodicityY:\n"
              << "    Cannot use y-periodicity with polar coordinates.\n";
    return;
  }
  if (s < Small) {
    std::cerr << m_className << "::SetPeriodicityY:\n"
              << "    Periodic length must be greater than zero.\n";
    return;
  }
 
  m_periodic[1] = true;
  m_sy = s;
  UpdatePeriodicity();
}

SetPolarCoordinates()

void Garfield::ComponentAnalyticField::SetPolarCoordinates

(

)

Use polar coordinates.

Definition at line 1593 of file ComponentAnalyticField.cc.

                                                 {
  if (!m_polar) {
    std::cout << m_className << "::SetPolarCoordinates:\n    "
              << "Switching to polar coordinates; resetting the cell.\n";
    CellInit();
  }
  m_polar = true;
  // Set default phi period.
  m_pery = true;
  m_sy = TwoPi;
}

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

                                                                {
  if (std::abs(xmax - xmin) < Small || std::abs(ymax - ymin) < Small) {
    std::cerr << m_className << "::SetScanningArea:\n"
              << "    Zero range not permitted.\n";
    return;
  }
  m_scanRange = ScanningRange::User;
  m_xScanMin = std::min(xmin, xmax);
  m_xScanMax = std::max(xmin, xmax);
  m_yScanMin = std::min(ymin, ymax);
  m_yScanMax = std::max(ymin, ymax);
}

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

                                                                         {
  m_scanRange = ScanningRange::FirstOrder;
  if (scale > 0.) {
    m_scaleRange = scale;
  } else {
    std::cerr << m_className << "::SetScanningAreaFirstOrder:\n"
              << "    Scaling factor must be > 0.\n";
  }
}

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

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

                                                                    {
  if (nX < 2) {
    std::cerr << m_className << "::SetScanningGrid:\n"
              << "    Number of x-lines must be > 1.\n";
  } else {
    m_nScanX = nX;
  }
  if (nY < 2) {
    std::cerr << m_className << "::SetScanningGrid:\n"
              << "    Number of y-lines must be > 1.\n";
  } else {
    m_nScanY = nY;
  }
}

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,z	coordinates [cm].
wx,wy,wz	components of the weighting field [1/cm].
label	name of the electrode

Reimplemented from Garfield::Component.

Definition at line 60 of file ComponentAnalyticField.hh.

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

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,z	coordinates [cm].
label	name of the electrode.

Returns

weighting potential [dimensionless].

Reimplemented from Garfield::Component.

Definition at line 67 of file ComponentAnalyticField.hh.

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

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

iw	index of the wire
detailed	flag to request a detailed calculation of the profile or only a fast one.
csag	coordinate along the wire.
xsag	x components of the sag profile.
ysag	y components of the sag profile.
stretch	relative elongation.
print	flag to print the calculation results or not.

Definition at line 2136 of file ComponentAnalyticField.cc.

                      {
  if (!m_cellset && !Prepare()) {
    std::cerr << m_className << "::WireDisplacement: Cell not set up.\n";
    return false;
  }
  if (iw >= m_nWires) {
    std::cerr << m_className
              << "::WireDisplacement: Wire index out of range.\n";
    return false;
  }
  if (m_polar) {
    std::cerr << m_className 
              << "::WireDisplacement: Cannot handle polar cells.\n";
    return false;
  }
  const auto& wire = m_w[iw];
  // Save the original coordinates.
  const double x0 = wire.x;
  const double y0 = wire.y;
  // First-order approximation.
  if (!Setup()) {
    std::cerr << m_className << "::WireDisplacement:\n"
              << "    Charge calculation failed at central position.\n";
    return false;
  }
 
  double fx = 0., fy = 0.;
  if (m_useElectrostaticForce) {
    double ex = 0., ey = 0.;
    ElectricFieldAtWire(iw, ex, ey);
    fx -= ex * wire.e * Internal2Newton;
    fy -= ey * wire.e * Internal2Newton;
  }
  if (m_useGravitationalForce) {
    // Mass per unit length [kg / cm].
    const double m = 1.e-3 * wire.density * Pi * wire.r * wire.r;
    fx -= m_down[0] * GravitationalAcceleration * m;
    fy -= m_down[1] * GravitationalAcceleration * m;
  }
  const double u2 = wire.u * wire.u;
  const double shiftx =
      -125 * fx * u2 / (GravitationalAcceleration * wire.tension);
  const double shifty =
      -125 * fy * u2 / (GravitationalAcceleration * wire.tension);
  // Get the elongation from this.
  const double s = 4 * sqrt(shiftx * shiftx + shifty * shifty) / wire.u;
  double length = wire.u;
  if (s > Small) {
    const double t = sqrt(1 + s * s);
    length *= 0.5 * (t + log(s + t) / s);
  }
  stretch = (length - wire.u) / wire.u;
  if (print) {
    std::cout << m_className
              << "::WireDisplacement:\n"
              << "    Forces and displacement in 0th order.\n"
              << "    Wire information: number   = " << iw << "\n"
              << "                      type     = " << wire.type << "\n"
              << "                      location = (" << wire.x << ", " << wire.y
              << ") cm\n"
              << "                      voltage  = " << wire.v << " V\n"
              << "                      length   = " << wire.u << " cm\n"
              << "                      tension  = " << wire.tension << " g\n"
              << "    In this position: Fx       = " << fx << " N/cm\n"
              << "                      Fy       = " << fy << " N/cm\n"
              << "                      x-shift  = " << shiftx << " cm\n"
              << "                      y-shift  = " << shifty << " cm\n"
              << "                      stretch  = " << 100. * stretch << "%\n";
  }
  if (!detailed) {
    csag = {0.};
    xsag = {shiftx};
    ysag = {shifty};
    return true;
  }
  std::vector<double> xMap(m_nScanX, 0.);
  std::vector<double> yMap(m_nScanY, 0.);
  std::vector<std::vector<double> > fxMap(m_nScanX,
                                          std::vector<double>(m_nScanY, 0.));
  std::vector<std::vector<double> > fyMap(m_nScanX,
                                          std::vector<double>(m_nScanY, 0.));
  if (!ForcesOnWire(iw, xMap, yMap, fxMap, fyMap)) {
    std::cerr << m_className << "::WireDisplacement:\n"
              << "    Could not compute the electrostatic force table.\n";
    return false;
  }
  // Compute the detailed wire shift.
  if (!SagDetailed(wire, xMap, yMap, fxMap, fyMap, csag, xsag, ysag)) {
    std::cerr << m_className << "::WireDisplacement: Wire " << iw << ".\n"
              << "    Computation of the wire sag failed.\n";
    return false;
  }
  // Verify that the wire is in range.
  const double sxmin = xMap.front();
  const double sxmax = xMap.back();
  const double symin = yMap.front();
  const double symax = yMap.back();
  const unsigned int nSag = xsag.size();
  bool outside = false;
  length = 0.;
  double xAvg = 0.;
  double yAvg = 0.;
  double xMax = 0.;
  double yMax = 0.;
  for (unsigned int i = 0; i < nSag; ++i) {
    if (x0 + xsag[i] < sxmin || x0 + xsag[i] > sxmax || y0 + ysag[i] < symin ||
        y0 + ysag[i] > symax) {
      outside = true;
    }
    xAvg += xsag[i];
    yAvg += ysag[i];
    xMax = std::max(xMax, std::abs(xsag[i]));
    yMax = std::max(yMax, std::abs(ysag[i]));
    if (i == 0) continue;
    const double dx = xsag[i] - xsag[i - 1];
    const double dy = ysag[i] - ysag[i - 1];
    const double dz = csag[i] - csag[i - 1];
    length += sqrt(dx * dx + dy * dy + dz * dz);
  }
  xAvg /= nSag;
  yAvg /= nSag;
  stretch = (length - wire.u) / wire.u;
  // Warn if a point outside the scanning area was found.
  if (outside) {
    std::cerr
        << m_className << "::WireDisplacement: Warning.\n    "
        << "The wire profile is located partially outside the scanning area.\n";
  }
  if (print) {
    std::cout << "    Sag profile for wire " << iw << ".\n"
              << " Point     z [cm]   x-sag [um]   y-sag [um]\n";
    for (unsigned int i = 0; i < nSag; ++i) {
      std::printf(" %3d   %10.4f  %10.4f  %10.4f\n", 
                  i, csag[i], xsag[i] * 1.e4, ysag[i] * 1.e4);
    }
    std::printf("    Average sag in x and y: %10.4f and %10.4f micron\n",
                1.e4 * xAvg, 1.e4 * yAvg);
    std::printf("    Maximum sag in x and y: %10.4f and %10.4f micron\n",
                1.e4 * xMax, 1.e4 * yMax);
    std::cout << "    Elongation:             " << 100. * stretch << "%\n";
  }
  return true;
}

---

The documentation for this class was generated from the following files:

• ComponentAnalyticField.hh
• ComponentAnalyticField.cc