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

Interpolation in a three-dimensional field map created by Sentaurus Device. More...

#include <ComponentTcad3d.hh>

+ Inheritance diagram for Garfield::ComponentTcad3d:

Public Member Functions

 ComponentTcad3d ()
 Constructor.
 
 ~ComponentTcad3d ()
 Destructor.
 
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).
 
void ElectricField (const double x, const double y, const double z, double &ex, double &ey, double &ez, Medium *&m, int &status) override
 
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
 
Medium * GetMedium (const double x, const double y, const double z) override
 Get the medium at a given location (x, y, z).
 
bool GetVoltageRange (double &vmin, double &vmax) override
 Calculate the voltage range [V].
 
bool GetBoundingBox (double &xmin, double &ymin, double &zmin, double &xmax, double &ymax, double &zmax) override
 Get the bounding box coordinates.
 
bool Initialise (const std::string &gridfilename, const std::string &datafilename)
 
bool SetWeightingField (const std::string &datfile1, const std::string &datfile2, const double dv)
 
void PrintRegions ()
 List all currently defined regions.
 
unsigned int GetNumberOfRegions () const
 Get the number of regions in the device.
 
void GetRegion (const unsigned int ireg, std::string &name, bool &active) const
 
void SetDriftRegion (const unsigned int ireg)
 
void UnsetDriftRegion (const unsigned int ireg)
 
void SetMedium (const unsigned int ireg, Medium *m)
 Set the medium for a given region.
 
bool GetMedium (const unsigned int ireg, Medium *&m) const
 Get the medium for a given region.
 
int GetNumberOfElements () const
 
bool GetElement (const unsigned int i, double &vol, double &dmin, double &dmax, int &type) const
 
bool GetElement (const unsigned int i, double &vol, double &dmin, double &dmax, int &type, int &node1, int &node2, int &node3, int &node4, int &node5, int &node6, int &node7, int &reg) const
 
unsigned int GetNumberOfNodes () const
 
bool GetNode (const unsigned int i, double &x, double &y, double &z, double &v, double &ex, double &ey, double &ez) const
 
- Public Member Functions inherited from Garfield::ComponentBase
 ComponentBase ()
 Constructor.
 
virtual ~ComponentBase ()
 Destructor.
 
virtual void SetGeometry (GeometryBase *geo)
 Define the geometry.
 
virtual void Clear ()
 Reset.
 
virtual Medium * GetMedium (const double x, const double y, const double z)
 Get the medium at a given location (x, y, z).
 
virtual void ElectricField (const double x, const double y, const double z, double &ex, double &ey, double &ez, Medium *&m, int &status)=0
 
virtual void ElectricField (const double x, const double y, const double z, double &ex, double &ey, double &ez, double &v, Medium *&m, int &status)=0
 Calculate the drift field [V/cm] and potential [V] at (x, y, z).
 
virtual bool GetVoltageRange (double &vmin, double &vmax)=0
 Calculate the voltage range [V].
 
virtual void WeightingField (const double x, const double y, const double z, double &wx, double &wy, double &wz, const std::string &label)
 
virtual double WeightingPotential (const double x, const double y, const double z, const std::string &label)
 
virtual void DelayedWeightingField (const double x, const double y, const double z, const double t, double &wx, double &wy, double &wz, const std::string &label)
 
virtual void MagneticField (const double x, const double y, const double z, double &bx, double &by, double &bz, int &status)
 
void SetMagneticField (const double bx, const double by, const double bz)
 Set a constant magnetic field.
 
virtual bool IsReady ()
 Ready for use?
 
virtual bool GetBoundingBox (double &xmin, double &ymin, double &zmin, double &xmax, double &ymax, double &zmax)
 Get the bounding box coordinates.
 
double IntegrateFluxCircle (const double xc, const double yc, const double r, const unsigned int nI=50)
 
double IntegrateFluxSphere (const double xc, const double yc, const double zc, const double r, const unsigned int nI=20)
 
double IntegrateFlux (const double x0, const double y0, const double z0, const double dx1, const double dy1, const double dz1, const double dx2, const double dy2, const double dz2, const unsigned int nU=20, const unsigned int nV=20)
 
virtual bool IsWireCrossed (const double x0, const double y0, const double z0, const double x1, const double y1, const double z1, double &xc, double &yc, double &zc, const bool centre, double &rc)
 
virtual bool IsInTrapRadius (const double q0, const double x0, const double y0, const double z0, double &xw, double &yw, double &rw)
 
void EnablePeriodicityX (const bool on=true)
 Enable simple periodicity in the $x$ direction.
 
void DisablePeriodicityX ()
 
void EnablePeriodicityY (const bool on=true)
 Enable simple periodicity in the $y$ direction.
 
void DisablePeriodicityY ()
 
void EnablePeriodicityZ (const bool on=true)
 Enable simple periodicity in the $z$ direction.
 
void DisablePeriodicityZ ()
 
void EnableMirrorPeriodicityX (const bool on=true)
 Enable mirror periodicity in the $x$ direction.
 
void DisableMirrorPeriodicityX ()
 
void EnableMirrorPeriodicityY (const bool on=true)
 Enable mirror periodicity in the $y$ direction.
 
void DisableMirrorPeriodicityY ()
 
void EnableMirrorPeriodicityZ (const bool on=true)
 Enable mirror periodicity in the $y$ direction.
 
void DisableMirrorPeriodicityZ ()
 
void EnableAxialPeriodicityX (const bool on=true)
 Enable axial periodicity in the $x$ direction.
 
void DisableAxialPeriodicityX ()
 
void EnableAxialPeriodicityY (const bool on=true)
 Enable axial periodicity in the $y$ direction.
 
void DisableAxialPeriodicityY ()
 
void EnableAxialPeriodicityZ (const bool on=true)
 Enable axial periodicity in the $z$ direction.
 
void DisableAxialPeriodicityZ ()
 
void EnableRotationSymmetryX (const bool on=true)
 Enable rotation symmetry around the $x$ axis.
 
void DisableRotationSymmetryX ()
 
void EnableRotationSymmetryY (const bool on=true)
 Enable rotation symmetry around the $y$ axis.
 
void DisableRotationSymmetryY ()
 
void EnableRotationSymmetryZ (const bool on=true)
 Enable rotation symmetry around the $z$ axis.
 
void DisableRotationSymmetryZ ()
 
void EnableDebugging ()
 Switch on debugging messages.
 
void DisableDebugging ()
 Switch off debugging messages.
 
void ActivateTraps ()
 Request trapping to be taken care of by the component (for TCAD).
 
void DeactivateTraps ()
 
bool IsTrapActive ()
 
void ActivateVelocityMap ()
 Request velocity to be taken care of by the component (for TCAD).
 
void DectivateVelocityMap ()
 
bool IsVelocityActive ()
 
virtual bool ElectronAttachment (const double, const double, const double, double &eta)
 Get the electron attachment coefficient.
 
virtual bool HoleAttachment (const double, const double, const double, double &eta)
 Get the hole attachment coefficient.
 
virtual void ElectronVelocity (const double, const double, const double, double &vx, double &vy, double &vz, Medium *&, int &status)
 Get the electron drift velocity.
 
virtual void HoleVelocity (const double, const double, const double, double &vx, double &vy, double &vz, Medium *&, int &status)
 Get the hole drift velocity.
 
virtual bool GetElectronLifetime (const double, const double, const double, double &etau)
 
virtual bool GetHoleLifetime (const double, const double, const double, double &htau)
 

Additional Inherited Members

virtual void Reset ()=0
 Reset the component.
 
virtual void UpdatePeriodicity ()=0
 Verify periodicities.
 
- Protected Attributes inherited from Garfield::ComponentBase
std::string m_className = "ComponentBase"
 Class name.
 
GeometryBase * m_geometry = nullptr
 Pointer to the geometry.
 
bool m_ready = false
 Ready for use?
 
bool m_activeTraps = false
 Does the component have traps?
 
bool m_hasVelocityMap = false
 Does the component have velocity maps?
 
std::array< bool, 3 > m_periodic = {{false, false, false}}
 Simple periodicity in x, y, z.
 
std::array< bool, 3 > m_mirrorPeriodic = {{false, false, false}}
 Mirror periodicity in x, y, z.
 
std::array< bool, 3 > m_axiallyPeriodic = {{false, false, false}}
 Axial periodicity in x, y, z.
 
std::array< bool, 3 > m_rotationSymmetric = {{false, false, false}}
 Rotation symmetry around x-axis, y-axis, z-axis.
 
double m_bx0 = 0.
 
double m_by0 = 0.
 
double m_bz0 = 0.
 
bool m_debug = false
 Switch on/off debugging messages.
 

Detailed Description

Interpolation in a three-dimensional field map created by Sentaurus Device.

Definition at line 12 of file ComponentTcad3d.hh.

Constructor & Destructor Documentation

◆ ComponentTcad3d()

Garfield::ComponentTcad3d::ComponentTcad3d ( )

Constructor.

Definition at line 14 of file ComponentTcad3d.cc.

14 : ComponentBase() {
15 m_className = "ComponentTcad3d";
16
17 m_regions.reserve(10);
18 m_vertices.reserve(10000);
19 m_elements.reserve(10000);
20}
Garfield::ComponentBase::ComponentBase
ComponentBase()
Constructor.
Definition: ComponentBase.cc:9
Garfield::ComponentBase::m_className
std::string m_className
Class name.
Definition: ComponentBase.hh:285

◆ ~ComponentTcad3d()

Garfield::ComponentTcad3d::~ComponentTcad3d ( )
inline

Destructor.

Definition at line 17 of file ComponentTcad3d.hh.

17{}

Member Function Documentation

◆ ElectricField() [1/2]

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

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

Implements Garfield::ComponentBase.

Definition at line 22 of file ComponentTcad3d.cc.

25 {
26 ex = ey = ez = p = 0.;
27 m = nullptr;
28 // Make sure the field map has been loaded.
29 if (!m_ready) {
30 std::cerr << m_className << "::ElectricField:\n"
31 << " Field map is not available for interpolation.\n";
32 status = -10;
33 return;
34 }
35
36 double x = xin, y = yin, z = zin;
37 // In case of periodicity, reduce to the cell volume.
38 bool xmirr = false, ymirr = false, zmirr = false;
39 MapCoordinates(x, y, z, xmirr, ymirr, zmirr);
40
41 // Check if the point is inside the bounding box.
42 if (x < m_xMinBB || x > m_xMaxBB || y < m_yMinBB || y > m_yMaxBB ||
43 z < m_zMinBB || z > m_zMaxBB) {
44 status = -6;
45 return;
46 }
47
48 // Assume this will work.
49 status = 0;
50 std::array<double, nMaxVertices> w;
51 const unsigned int i = FindElement(x, y, z, w);
52 if (i >= m_elements.size()) {
53 // Point is outside the mesh.
54 status = -6;
55 return;
56 }
57 const Element& element = m_elements[i];
58 const unsigned int nVertices = element.type == 2 ? 3 : 4;
59 for (unsigned int j = 0; j < nVertices; ++j) {
60 const Vertex& vj = m_vertices[element.vertex[j]];
61 ex += w[j] * vj.ex;
62 ey += w[j] * vj.ey;
63 ez += w[j] * vj.ez;
64 p += w[j] * vj.p;
65 }
66 if (xmirr) ex = -ex;
67 if (ymirr) ey = -ey;
68 if (zmirr) ez = -ez;
69 m = m_regions[element.region].medium;
70 if (!m_regions[element.region].drift || !m) status = -5;
71 m_lastElement = i;
72}
Garfield::ComponentBase::m_ready
bool m_ready
Ready for use?
Definition: ComponentBase.hh:291

Referenced by ElectricField().

◆ ElectricField() [2/2]

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

Calculate the drift field at given point.

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

Status flags: 

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

Implements Garfield::ComponentBase.

Definition at line 74 of file ComponentTcad3d.cc.

76 {
77 double v = 0.;
78 ElectricField(x, y, z, ex, ey, ez, v, m, status);
79}
Garfield::ComponentTcad3d::ElectricField
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).
Definition: ComponentTcad3d.cc:22

◆ GetBoundingBox()

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

Get the bounding box coordinates.

Reimplemented from Garfield::ComponentBase.

Definition at line 483 of file ComponentTcad3d.cc.

484 {
485 if (!m_ready) return false;
486 xmin = m_xMinBB;
487 ymin = m_yMinBB;
488 zmin = m_zMinBB;
489 xmax = m_xMaxBB;
490 ymax = m_yMaxBB;
491 zmax = m_zMaxBB;
492 if (m_periodic[0] || m_mirrorPeriodic[0]) {
493 xmin = -INFINITY;
494 xmax = +INFINITY;
495 }
496 if (m_periodic[1] || m_mirrorPeriodic[1]) {
497 ymin = -INFINITY;
498 ymax = +INFINITY;
499 }
500 if (m_periodic[2] || m_mirrorPeriodic[2]) {
501 zmin = -INFINITY;
502 zmax = +INFINITY;
503 }
504 return true;
505}
Garfield::ComponentBase::m_periodic
std::array< bool, 3 > m_periodic
Simple periodicity in x, y, z.
Definition: ComponentBase.hh:299
Garfield::ComponentBase::m_mirrorPeriodic
std::array< bool, 3 > m_mirrorPeriodic
Mirror periodicity in x, y, z.
Definition: ComponentBase.hh:301

◆ GetElement() [1/2]

bool Garfield::ComponentTcad3d::GetElement ( const unsigned int  i,
double &  vol,
double &  dmin,
double &  dmax,
int &  type 
) const

Definition at line 641 of file ComponentTcad3d.cc.

642 {
643 if (i >= m_elements.size()) {
644 std::cerr << m_className << "::GetElement: Index out of range.\n";
645 return false;
646 }
647
648 const Element& element = m_elements[i];
649 if (element.type == 2) {
650 // Triangle
651 const Vertex& v0 = m_vertices[element.vertex[0]];
652 const Vertex& v1 = m_vertices[element.vertex[1]];
653 const Vertex& v2 = m_vertices[element.vertex[2]];
654 const double vx =
655 (v1.y - v0.y) * (v2.z - v0.z) - (v1.z - v0.z) * (v2.y - v0.y);
656 const double vy =
657 (v1.z - v0.z) * (v2.x - v0.x) - (v1.x - v0.x) * (v2.z - v0.z);
658 const double vz =
659 (v1.x - v0.x) * (v2.y - v0.y) - (v1.y - v0.y) * (v2.x - v0.x);
660 vol = sqrt(vx * vx + vy * vy + vz * vz);
661 const double a =
662 sqrt(pow(v1.x - v0.x, 2) + pow(v1.y - v0.y, 2) + pow(v1.z - v0.z, 2));
663 const double b =
664 sqrt(pow(v2.x - v0.x, 2) + pow(v2.y - v0.y, 2) + pow(v2.z - v0.z, 2));
665 const double c =
666 sqrt(pow(v1.x - v2.x, 2) + pow(v1.y - v2.y, 2) + pow(v1.z - v2.z, 2));
667 dmin = dmax = a;
668 if (b < dmin) dmin = b;
669 if (c < dmin) dmin = c;
670 if (b > dmax) dmax = b;
671 if (c > dmax) dmax = c;
672 } else if (element.type == 5) {
673 // Tetrahedron
674 const Vertex& v0 = m_vertices[element.vertex[0]];
675 const Vertex& v1 = m_vertices[element.vertex[1]];
676 const Vertex& v2 = m_vertices[element.vertex[2]];
677 const Vertex& v3 = m_vertices[element.vertex[3]];
678 vol = fabs((v3.x - v0.x) * ((v1.y - v0.y) * (v2.z - v0.z) -
679 (v2.y - v0.y) * (v1.z - v0.z)) +
680 (v3.y - v0.y) * ((v1.z - v0.z) * (v2.x - v0.x) -
681 (v2.z - v0.z) * (v1.x - v0.x)) +
682 (v3.z - v0.z) * ((v1.x - v0.x) * (v2.y - v0.y) -
683 (v3.x - v0.x) * (v1.y - v0.y))) /
684 6.;
685 // Loop over all pairs of m_vertices.
686 for (int j = 0; j < nMaxVertices - 1; ++j) {
687 const Vertex& vj = m_vertices[element.vertex[j]];
688 for (int k = j + 1; k < nMaxVertices; ++k) {
689 const Vertex& vk = m_vertices[element.vertex[k]];
690 // Compute distance.
691 const double dist = sqrt(pow(vj.x - vk.x, 2) + pow(vj.y - vk.y, 2) +
692 pow(vj.z - vk.z, 2));
693 if (k == 1) {
694 dmin = dmax = dist;
695 } else {
696 if (dist < dmin) dmin = dist;
697 if (dist > dmax) dmax = dist;
698 }
699 }
700 }
701 } else {
702 std::cerr << m_className << "::GetElement:\n"
703 << " Unexpected element type (" << type << ").\n";
704 return false;
705 }
706 return true;
707}
Heed::pow
DoubleAc pow(const DoubleAc &f, double p)
Definition: DoubleAc.cpp:337
Heed::fabs
DoubleAc fabs(const DoubleAc &f)
Definition: DoubleAc.h:615
Heed::sqrt
DoubleAc sqrt(const DoubleAc &f)
Definition: DoubleAc.cpp:314

Referenced by GetElement().

◆ GetElement() [2/2]

bool Garfield::ComponentTcad3d::GetElement ( const unsigned int  i,
double &  vol,
double &  dmin,
double &  dmax,
int &  type,
int &  node1,
int &  node2,
int &  node3,
int &  node4,
int &  node5,
int &  node6,
int &  node7,
int &  reg 
) const

Definition at line 709 of file ComponentTcad3d.cc.

713 {
714 if (!GetElement(i, vol, dmin, dmax, type)) return false;
715 const Element& element = m_elements[i];
716 node1 = element.vertex[0];
717 node2 = element.vertex[1];
718 node3 = element.vertex[2];
719 node4 = element.vertex[3];
720 node5 = element.vertex[4];
721 node6 = element.vertex[5];
722 node7 = element.vertex[6];
723 reg = element.region;
724 return true;
725}
Garfield::ComponentTcad3d::GetElement
bool GetElement(const unsigned int i, double &vol, double &dmin, double &dmax, int &type) const
Definition: ComponentTcad3d.cc:641

◆ GetMedium() [1/2]

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

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

Reimplemented from Garfield::ComponentBase.

Definition at line 145 of file ComponentTcad3d.cc.

146 {
147 // Make sure the field map has been loaded.
148 if (!m_ready) {
149 std::cerr << m_className << "::GetMedium:\n"
150 << " Field map not available for interpolation.\n";
151 return nullptr;
152 }
153
154 double x = xin, y = yin, z = zin;
155 bool xmirr = false, ymirr = false, zmirr = false;
156 MapCoordinates(x, y, z, xmirr, ymirr, zmirr);
157
158 // Check if the point is inside the bounding box.
159 if (x < m_xMinBB || x > m_xMaxBB || y < m_yMinBB || y > m_yMaxBB ||
160 z < m_zMinBB || z > m_zMaxBB) {
161 return nullptr;
162 }
163
164 std::array<double, nMaxVertices> w;
165 const unsigned int i = FindElement(x, y, z, w);
166 if (i >= m_elements.size()) {
167 // Point is outside the mesh.
168 return nullptr;
169 }
170 m_lastElement = i;
171 const Element& element = m_elements[i];
172 return m_regions[element.region].medium;
173}

◆ GetMedium() [2/2]

bool Garfield::ComponentTcad3d::GetMedium ( const unsigned int  ireg,
Medium *&  m 
) const

Get the medium for a given region.

Definition at line 586 of file ComponentTcad3d.cc.

586 {
587 if (i >= m_regions.size()) {
588 std::cerr << m_className << "::GetMedium: Index out of range.\n";
589 return false;
590 }
591
592 m = m_regions[i].medium;
593 if (!m) return false;
594 return true;
595}

◆ GetNode()

bool Garfield::ComponentTcad3d::GetNode ( const unsigned int  i,
double &  x,
double &  y,
double &  z,
double &  v,
double &  ex,
double &  ey,
double &  ez 
) const

Definition at line 727 of file ComponentTcad3d.cc.

729 {
730 if (i >= m_vertices.size()) {
731 std::cerr << m_className << "::GetNode: Index out of range.\n";
732 return false;
733 }
734
735 const Vertex& vi = m_vertices[i];
736 x = vi.x;
737 y = vi.y;
738 z = vi.z;
739 v = vi.p;
740 ex = vi.ex;
741 ey = vi.ey;
742 ez = vi.ez;
743 return true;
744}

◆ GetNumberOfElements()

int Garfield::ComponentTcad3d::GetNumberOfElements ( ) const
inline

Definition at line 69 of file ComponentTcad3d.hh.

69{ return m_elements.size(); }

◆ GetNumberOfNodes()

unsigned int Garfield::ComponentTcad3d::GetNumberOfNodes ( ) const
inline

Definition at line 75 of file ComponentTcad3d.hh.

75{ return m_vertices.size(); }

◆ GetNumberOfRegions()

unsigned int Garfield::ComponentTcad3d::GetNumberOfRegions ( ) const
inline

Get the number of regions in the device.

Definition at line 59 of file ComponentTcad3d.hh.

59{ return m_regions.size(); }

◆ GetRegion()

void Garfield::ComponentTcad3d::GetRegion ( const unsigned int  ireg,
std::string &  name,
bool &  active 
) const

Definition at line 547 of file ComponentTcad3d.cc.

548 {
549 if (i >= m_regions.size()) {
550 std::cerr << m_className << "::GetRegion: Index out of range.\n";
551 return;
552 }
553 name = m_regions[i].name;
554 active = m_regions[i].drift;
555}

◆ GetVoltageRange()

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

Calculate the voltage range [V].

Implements Garfield::ComponentBase.

Definition at line 507 of file ComponentTcad3d.cc.

507 {
508 if (!m_ready) return false;
509 vmin = m_pMin;
510 vmax = m_pMax;
511 return true;
512}

◆ Initialise()

bool Garfield::ComponentTcad3d::Initialise ( const std::string &  gridfilename,
const std::string &  datafilename 
)

Import mesh and field map from files.

Parameters
gridfilenamename of the .grd file containing the mesh
datafilenamename of the .dat file containing the nodal solution

Definition at line 175 of file ComponentTcad3d.cc.

176 {
177 m_ready = false;
178 // Import mesh data from .grd file.
179 if (!LoadGrid(gridfilename)) {
180 std::cerr << m_className << "::Initialise:\n"
181 << " Importing mesh data failed.\n";
182 return false;
183 }
184
185 // Import electric field and potential from .dat file.
186 if (!LoadData(datafilename)) {
187 std::cerr << m_className << "::Initialise:\n"
188 << " Importing electric field and potential failed.\n";
189 return false;
190 }
191
192 // Find min./max. coordinates and potentials.
193 m_xMaxBB = m_vertices[m_elements[0].vertex[0]].x;
194 m_yMaxBB = m_vertices[m_elements[0].vertex[0]].y;
195 m_zMaxBB = m_vertices[m_elements[0].vertex[0]].z;
196 m_xMinBB = m_xMaxBB;
197 m_yMinBB = m_yMaxBB;
198 m_zMinBB = m_zMaxBB;
199 m_pMax = m_pMin = m_vertices[m_elements[0].vertex[0]].p;
200 const unsigned int nElements = m_elements.size();
201 for (unsigned int i = 0; i < nElements; ++i) {
202 Element& element = m_elements[i];
203 double xmin = m_vertices[element.vertex[0]].x;
204 double ymin = m_vertices[element.vertex[0]].y;
205 double zmin = m_vertices[element.vertex[0]].z;
206 double xmax = xmin;
207 double ymax = ymin;
208 double zmax = zmin;
209 const unsigned int nVertices = element.type == 2 ? 3 : 4;
210 for (unsigned int j = 0; j < nVertices; ++j) {
211 const Vertex& vj = m_vertices[element.vertex[j]];
212 if (vj.x < xmin) xmin = vj.x;
213 if (vj.x > xmax) xmax = vj.x;
214 if (vj.y < ymin) ymin = vj.y;
215 if (vj.y > ymax) ymax = vj.y;
216 if (vj.z < zmin) zmin = vj.z;
217 if (vj.z > zmax) zmax = vj.z;
218 if (vj.p < m_pMin) m_pMin = vj.p;
219 if (vj.p > m_pMax) m_pMax = vj.p;
220 }
221 constexpr double tol = 1.e-6;
222 element.xmin = xmin - tol;
223 element.xmax = xmax + tol;
224 element.ymin = ymin - tol;
225 element.ymax = ymax + tol;
226 element.zmin = zmin - tol;
227 element.zmax = zmax + tol;
228 m_xMinBB = std::min(m_xMinBB, xmin);
229 m_xMaxBB = std::max(m_xMaxBB, xmax);
230 m_yMinBB = std::min(m_yMinBB, ymin);
231 m_yMaxBB = std::max(m_yMaxBB, ymax);
232 m_zMinBB = std::min(m_zMinBB, zmin);
233 m_zMaxBB = std::max(m_zMaxBB, zmax);
234 }
235
236 std::cout << m_className << "::Initialise:\n"
237 << " Bounding box:\n"
238 << " " << m_xMinBB << " < x [cm] < " << m_xMaxBB << "\n"
239 << " " << m_yMinBB << " < y [cm] < " << m_yMaxBB << "\n"
240 << " " << m_zMinBB << " < z [cm] < " << m_zMaxBB << "\n"
241 << " Voltage range:\n"
242 << " " << m_pMin << " < V < " << m_pMax << "\n";
243
244 bool ok = true;
245
246 // Count the number of elements belonging to a region.
247 const int nRegions = m_regions.size();
248 std::vector<unsigned int> nElementsRegion(nRegions, 0);
249
250 // Count the different element shapes.
251 unsigned int nTriangles = 0;
252 unsigned int nTetrahedra = 0;
253 unsigned int nOtherShapes = 0;
254
255 // Check if there are elements which are not part of any region.
256 unsigned int nLoose = 0;
257 std::vector<int> looseElements;
258
259 // Check if there are degenerate elements.
260 unsigned int nDegenerate = 0;
261 std::vector<int> degenerateElements;
262
263 for (unsigned int i = 0; i < nElements; ++i) {
264 const Element& element = m_elements[i];
265 if (element.type == 2) {
266 ++nTriangles;
267 if (element.vertex[0] == element.vertex[1] ||
268 element.vertex[1] == element.vertex[2] ||
269 element.vertex[2] == element.vertex[0]) {
270 degenerateElements.push_back(i);
271 ++nDegenerate;
272 }
273 } else if (element.type == 5) {
274 if (element.vertex[0] == element.vertex[1] ||
275 element.vertex[0] == element.vertex[2] ||
276 element.vertex[0] == element.vertex[3] ||
277 element.vertex[1] == element.vertex[2] ||
278 element.vertex[1] == element.vertex[3] ||
279 element.vertex[2] == element.vertex[3]) {
280 degenerateElements.push_back(i);
281 ++nDegenerate;
282 }
283 ++nTetrahedra;
284 } else {
285 // Other shapes should not occur, since they were excluded in LoadGrid.
286 ++nOtherShapes;
287 }
288 if (element.region >= 0 && element.region < nRegions) {
289 ++nElementsRegion[element.region];
290 } else {
291 looseElements.push_back(i);
292 ++nLoose;
293 }
294 }
295
296 if (nDegenerate > 0) {
297 std::cerr << m_className << "::Initialise:\n"
298 << " The following elements are degenerate:\n";
299 for (unsigned int i = 0; i < nDegenerate; ++i) {
300 std::cerr << " " << degenerateElements[i] << "\n";
301 }
302 ok = false;
303 }
304
305 if (nLoose > 0) {
306 std::cerr << m_className << "::Initialise:\n"
307 << " The following elements are not part of any region:\n";
308 for (unsigned int i = 0; i < nLoose; ++i) {
309 std::cerr << " " << looseElements[i] << "\n";
310 }
311 ok = false;
312 }
313
314 std::cout << m_className << "::Initialise:\n"
315 << " Number of regions: " << nRegions << "\n";
316 for (int i = 0; i < nRegions; ++i) {
317 std::cout << " " << i << ": " << m_regions[i].name << ", "
318 << nElementsRegion[i] << " elements\n";
319 }
320
321 std::cout << " Number of elements: " << nElements << "\n";
322 if (nTriangles > 0) {
323 std::cout << " " << nTriangles << " triangles\n";
324 }
325 if (nTetrahedra > 0) {
326 std::cout << " " << nTetrahedra << " tetrahedra\n";
327 }
328 if (nOtherShapes > 0) {
329 std::cerr << " " << nOtherShapes << " elements of unknown type\n"
330 << " Program bug!\n";
331 m_ready = false;
332 Cleanup();
333 return false;
334 }
335 if (m_debug) {
336 // For each element, print the indices of the constituting vertices.
337 for (unsigned int i = 0; i < nElements; ++i) {
338 const Element& element = m_elements[i];
339 if (element.type == 2) {
340 std::cout << " " << i << ": " << element.vertex[0] << " "
341 << element.vertex[1] << " " << element.vertex[2]
342 << " (triangle, region " << element.region << ")\n";
343 } else if (element.type == 5) {
344 std::cout << " " << i << ": " << element.vertex[0] << " "
345 << element.vertex[1] << " " << element.vertex[2] << " "
346 << element.vertex[3] << " (tetrahedron, region "
347 << element.region << ")\n";
348 }
349 }
350 }
351
352 const unsigned int nVertices = m_vertices.size();
353 std::cout << " Number of vertices: " << nVertices << "\n";
354 if (m_debug) {
355 for (unsigned int i = 0; i < nVertices; ++i) {
356 const Vertex& vi = m_vertices[i];
357 std::cout << " " << i << ": (x, y, z) = (" << vi.x << ", " << vi.y
358 << ", " << vi.z << "), V = " << vi.p << "\n";
359 }
360 }
361
362 // Find adjacent elements.
363 std::cout << m_className << "::Initialise:\n"
364 << " Looking for neighbouring elements. Be patient...\n";
365 FindNeighbours();
366
367 if (!ok) {
368 m_ready = false;
369 Cleanup();
370 return false;
371 }
372
373 m_ready = true;
374 UpdatePeriodicity();
375 std::cout << m_className << "::Initialise: Initialisation finished.\n";
376 return true;
377}
Garfield::ComponentBase::m_debug
bool m_debug
Switch on/off debugging messages.
Definition: ComponentBase.hh:310

◆ PrintRegions()

void Garfield::ComponentTcad3d::PrintRegions ( )

List all currently defined regions.

Definition at line 514 of file ComponentTcad3d.cc.

514 {
515 // Do not proceed if not properly initialised.
516 if (!m_ready) {
517 std::cerr << m_className << "::PrintRegions:\n"
518 << " Field map not yet initialised.\n";
519 return;
520 }
521
522 if (m_regions.empty()) {
523 std::cerr << m_className << "::PrintRegions:\n"
524 << " No regions are currently defined.\n";
525 return;
526 }
527
528 const unsigned int nRegions = m_regions.size();
529 std::cout << m_className << "::PrintRegions:\n"
530 << " Currently " << nRegions << " regions are defined.\n"
531 << " Index Name Medium\n";
532 for (unsigned int i = 0; i < nRegions; ++i) {
533 std::cout << " " << i << " " << m_regions[i].name;
534 if (!m_regions[i].medium) {
535 std::cout << " none ";
536 } else {
537 std::cout << " " << m_regions[i].medium->GetName();
538 }
539 if (m_regions[i].drift) {
540 std::cout << " (active region)\n";
541 } else {
542 std::cout << "\n";
543 }
544 }
545}

◆ SetDriftRegion()

void Garfield::ComponentTcad3d::SetDriftRegion ( const unsigned int  ireg)

Definition at line 557 of file ComponentTcad3d.cc.

557 {
558 if (i >= m_regions.size()) {
559 std::cerr << m_className << "::SetDriftRegion: Index out of range.\n";
560 return;
561 }
562 m_regions[i].drift = true;
563}

◆ SetMedium()

void Garfield::ComponentTcad3d::SetMedium ( const unsigned int  ireg,
Medium *  m 
)

Set the medium for a given region.

Definition at line 573 of file ComponentTcad3d.cc.

573 {
574 if (i >= m_regions.size()) {
575 std::cerr << m_className << "::SetMedium: Index out of range.\n";
576 return;
577 }
578
579 if (!medium) {
580 std::cerr << m_className << "::SetMedium: Null pointer.\n";
581 return;
582 }
583 m_regions[i].medium = medium;
584}

◆ SetWeightingField()

bool Garfield::ComponentTcad3d::SetWeightingField ( const std::string &  datfile1,
const std::string &  datfile2,
const double  dv 
)

Import field maps defining the weighting field and potential.

Parameters
datfile1.dat file containing the field map at nominal bias.
datfile2.dat file containing the field map for a configuration with the potential at the electrode to be read out increased by a small voltage dv.
dvincrease in electrode potential between the two field maps.

The field maps must use the same mesh as the drift field.

Definition at line 379 of file ComponentTcad3d.cc.

381 {
382
383 if (!m_ready) {
384 std::cerr << m_className << "::SetWeightingField:\n"
385 << " Mesh is not available. Call Initialise first.\n";
386 return false;
387 }
388 if (dv < Small) {
389 std::cerr << m_className << "::SetWeightingField:\n"
390 << " Voltage difference must be > 0.\n";
391 return false;
392 }
393 const double s = 1. / dv;
394
395 m_wf.clear();
396 m_wp.clear();
397 // Load first the field/potential at nominal bias.
398 std::vector<std::array<double, 3> > wf1;
399 std::vector<double> wp1;
400 if (!LoadWeightingField(datfile1, wf1, wp1)) {
401 std::cerr << m_className << "::SetWeightingField:\n"
402 << " Could not import data from " << datfile1 << ".\n";
403 return false;
404 }
405 // Then load the field/potential for the configuration with the potential
406 // at the electrode to be read out increased by small voltage dv.
407 std::vector<std::array<double, 3> > wf2;
408 std::vector<double> wp2;
409 if (!LoadWeightingField(datfile2, wf2, wp2)) {
410 std::cerr << m_className << "::SetWeightingField:\n"
411 << " Could not import data from " << datfile2 << ".\n";
412 return false;
413 }
414 const unsigned int nVertices = m_vertices.size();
415 bool foundField = true;
416 if (wf1.size() != nVertices || wf2.size() != nVertices) {
417 foundField = false;
418 std::cerr << m_className << "::SetWeightingField:\n"
419 << " Could not load electric field values.\n";
420 }
421 bool foundPotential = true;
422 if (wp1.size() != nVertices || wp2.size() != nVertices) {
423 foundPotential = false;
424 std::cerr << m_className << "::SetWeightingField:\n"
425 << " Could not load electrostatic potentials.\n";
426 }
427 if (!foundField && !foundPotential) return false;
428 if (foundField) {
429 m_wf.assign(nVertices, {0., 0., 0.});
430 for (unsigned int i = 0; i < nVertices; ++i) {
431 for (unsigned int j = 0; j < 3; ++j) {
432 m_wf[i][j] = (wf2[i][j] - wf1[i][j]) * s;
433 }
434 }
435 }
436 if (foundPotential) {
437 m_wp.assign(nVertices, 0.);
438 for (unsigned int i = 0; i < nVertices; ++i) {
439 m_wp[i] = (wp2[i] - wp1[i]) * s;
440 }
441 }
442 return true;
443}

◆ UnsetDriftRegion()

void Garfield::ComponentTcad3d::UnsetDriftRegion ( const unsigned int  ireg)

Definition at line 565 of file ComponentTcad3d.cc.

565 {
566 if (i >= m_regions.size()) {
567 std::cerr << m_className << "::UnsetDriftRegion: Index out of range.\n";
568 return;
569 }
570 m_regions[i].drift = false;
571}

◆ WeightingField()

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

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

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

Reimplemented from Garfield::ComponentBase.

Definition at line 81 of file ComponentTcad3d.cc.

83 {
84 wx = wy = wz = 0.;
85 if (m_wf.empty()) {
86 std::cerr << m_className << "::WeightingField: Not available.\n";
87 return;
88 }
89 double x = xin, y = yin, z = zin;
90 // In case of periodicity, reduce to the cell volume.
91 bool xmirr = false, ymirr = false, zmirr = false;
92 MapCoordinates(x, y, z, xmirr, ymirr, zmirr);
93
94 // Check if the point is inside the bounding box.
95 if (x < m_xMinBB || x > m_xMaxBB || y < m_yMinBB || y > m_yMaxBB ||
96 z < m_zMinBB || z > m_zMaxBB) return;
97
98 std::array<double, nMaxVertices> w;
99 const unsigned int i = FindElement(x, y, z, w);
100 if (i >= m_elements.size()) return;
101
102 const Element& element = m_elements[i];
103 const unsigned int nVertices = element.type == 2 ? 3 : 4;
104 for (unsigned int j = 0; j < nVertices; ++j) {
105 const auto& f = m_wf[element.vertex[j]];
106 wx += w[j] * f[0];
107 wy += w[j] * f[1];
108 wz += w[j] * f[2];
109 }
110 if (xmirr) wx = -wx;
111 if (ymirr) wy = -wy;
112 if (zmirr) wz = -wz;
113}

◆ WeightingPotential()

double Garfield::ComponentTcad3d::WeightingPotential ( const double  x,
const double  y,
const double  z,
const std::string &  label 
)
overridevirtual

Calculate the weighting potential at a given point.

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

Reimplemented from Garfield::ComponentBase.

Definition at line 115 of file ComponentTcad3d.cc.

117 {
118
119 if (m_wp.empty()) {
120 std::cerr << m_className << "::WeightingPotential: Not available.\n";
121 return 0.;
122 }
123 // In case of periodicity, reduce to the cell volume.
124 double x = xin, y = yin, z = zin;
125 bool xmirr = false, ymirr = false, zmirr = false;
126 MapCoordinates(x, y, z, xmirr, ymirr, zmirr);
127
128 // Check if the point is inside the bounding box.
129 if (x < m_xMinBB || x > m_xMaxBB || y < m_yMinBB || y > m_yMaxBB ||
130 z < m_zMinBB || z > m_zMaxBB) return 0.;
131
132 std::array<double, nMaxVertices> w;
133 const unsigned int i = FindElement(x, y, z, w);
134 if (i >= m_elements.size()) return 0.;
135
136 double v = 0.;
137 const Element& element = m_elements[i];
138 const unsigned int nVertices = element.type == 2 ? 3 : 4;
139 for (unsigned int j = 0; j < nVertices; ++j) {
140 v += w[j] * m_wp[element.vertex[j]];
141 }
142 return v;
143}
The documentation for this class was generated from the following files: