Garfield++ v2r0
A toolkit for the detailed simulation of particle detectors based on ionisation measurement in gases and semiconductors
Loading...
Searching...
No Matches
ComponentAnalyticField.hh
Go to the documentation of this file.
1#ifndef G_COMPONENT_ANALYTIC_FIELD_H
2#define G_COMPONENT_ANALYTIC_FIELD_H
3
4#include <cmath>
5#include <complex>
6
7#include "ComponentBase.hh"
9
10namespace Garfield {
11
12/// Semi-analytic calculation of two-dimensional configurations
13/// consisting of wires, planes, and tubes.
14
16
17 public:
18 /// Constructor
20 /// Destructor
22
23 void ElectricField(const double x, const double y, const double z, double& ex,
24 double& ey, double& ez, Medium*& m, int& status) {
25
26 m = NULL;
27 // Calculate the field.
28 double v = 0.;
29 status = Field(x, y, z, ex, ey, ez, v, false);
30 // If the field is ok, get the medium.
31 if (status == 0) {
32 m = GetMedium(x, y, z);
33 if (!m) {
34 status = -6;
35 } else if (!m->IsDriftable()) {
36 status = -5;
37 }
38 }
39 }
40
41 void ElectricField(const double x, const double y, const double z, double& ex,
42 double& ey, double& ez, double& v, Medium*& m,
43 int& status) {
44 m = NULL;
45 // Calculate the field.
46 status = Field(x, y, z, ex, ey, ez, v, true);
47 // If the field is ok, get the medium.
48 if (status == 0) {
49 m = GetMedium(x, y, z);
50 if (!m) {
51 status = -6;
52 } else if (!m->IsDriftable()) {
53 status = -5;
54 }
55 }
56 }
57
58 bool GetVoltageRange(double& pmin, double& pmax);
59
60 void WeightingField(const double x, const double y, const double z,
61 double& wx, double& wy, double& wz,
62 const std::string& label) {
63 wx = wy = wz = 0.;
64 double volt = 0.;
65 if (!m_sigset) PrepareSignals();
66 Wfield(x, y, z, wx, wy, wz, volt, label, false);
67 }
68 double WeightingPotential(const double x, const double y, const double z,
69 const std::string& label) {
70 double wx = 0., wy = 0., wz = 0.;
71 double volt = 0.;
72 if (!m_sigset) PrepareSignals();
73 Wfield(x, y, z, wx, wy, wz, volt, label, true);
74 return volt;
75 }
76
77 bool GetBoundingBox(double& x0, double& y0, double& z0, double& x1,
78 double& y1, double& z1);
79
80 bool IsWireCrossed(const double x0, const double y0, const double z0,
81 const double x1, const double y1, const double z1,
82 double& xc, double& yc, double& zc);
83
84 bool IsInTrapRadius(const double q0, const double x0, const double y0,
85 const double z0, double& xw, double& yx, double& rw);
86
87 /// Add a wire at (x, y) .
88 void AddWire(const double x, const double y, const double diameter,
89 const double voltage, const std::string& label,
90 const double length = 100., const double tension = 50.,
91 const double rho = 19.3, const int ntrap = 5);
92 /// Add a tube.
93 void AddTube(const double radius, const double voltage, const int nEdges,
94 const std::string& label);
95 /// Add a plane at constant x.
96 void AddPlaneX(const double x, const double voltage,
97 const std::string& label);
98 /// Add a plane at constant y.
99 void AddPlaneY(const double y, const double voltage,
100 const std::string& label);
101
102 void AddStripOnPlaneX(const char direction, const double x, const double smin,
103 const double smax, const std::string& label,
104 const double gap = -1.);
105 void AddStripOnPlaneY(const char direction, const double y, const double smin,
106 const double smax, const std::string& label,
107 const double gap = -1.);
108 void AddPixelOnPlaneX(const double x, const double ymin, const double ymax,
109 const double zmin, const double zmax,
110 const std::string& label, const double gap = -1.);
111 void AddPixelOnPlaneY(const double y, const double xmin, const double xmax,
112 const double zmin, const double zmax,
113 const std::string& label, const double gap = -1.);
114
115 /// Set the periodic length [cm] in the x-direction.
116 void SetPeriodicityX(const double s);
117 /// Set the periodic length [cm] in the y-direction.
118 void SetPeriodicityY(const double s);
119 bool GetPeriodicityX(double& s);
120 bool GetPeriodicityY(double& s);
121
122 /// Add a point charge.
123 void AddCharge(const double x, const double y, const double z,
124 const double q);
125 void ClearCharges();
126 void PrintCharges() const;
127
128 /** Return the cell type.
129 * Cells are classified according to the number
130 * and orientation of planes, the presence of
131 * periodicities and the location of the wires
132 * as one of the following types:
133 *
134 * A non-periodic cells with at most 1 x- and 1 y-plane
135 * B1X x-periodic cells without x-planes and at most 1 y-plane
136 * B1Y y-periodic cells without y-planes and at most 1 x-plane
137 * B2X cells with 2 x-planes and at most 1 y-plane
138 * B2Y cells with 2 y-planes and at most 1 x-plane
139 * C1 doubly periodic cells without planes
140 * C2X doubly periodic cells with x-planes
141 * C2Y doubly periodic cells with y-planes
142 * C3 double periodic cells with x- and y-planes
143 * D1 round tubes without axial periodicity
144 * D2 round tubes with axial periodicity
145 * D3 polygonal tubes without axial periodicity
146 */
147 std::string GetCellType() {
148 if (!m_cellset) {
149 if (CellCheck()) CellType();
150 }
151 return m_scellType;
152 }
153
154 /// Setup the weighting field for a given group of wires or planes.
155 void AddReadout(const std::string& label);
156
157 void EnableChargeCheck(const bool on = true) { m_chargeCheck = on; }
159
160 unsigned int GetNumberOfWires() const { return m_nWires; }
161 bool GetWire(const unsigned int i, double& x, double& y, double& diameter,
162 double& voltage, std::string& label, double& length,
163 double& charge, int& ntrap) const;
164
165 unsigned int GetNumberOfPlanesX() const;
166 unsigned int GetNumberOfPlanesY() const;
167 bool GetPlaneX(const unsigned int i, double& x, double& voltage,
168 std::string& label) const;
169 bool GetPlaneY(const unsigned int i, double& y, double& voltage,
170 std::string& label) const;
171
172 bool GetTube(double& r, double& voltage, int& nEdges, std::string& label) const;
173
174 enum Cell {
188 Unknown
189 };
190
191 private:
192 bool m_chargeCheck;
193
194 bool m_cellset;
195 bool m_sigset;
196
197 bool m_polar;
198
199 // Cell type (as string and number)
200 std::string m_scellType;
201 Cell m_cellType;
202
203 // Bounding box
204 double m_xmin, m_xmax;
205 double m_ymin, m_ymax;
206 double m_zmin, m_zmax;
207
208 // Voltage range
209 double vmin, vmax;
210
211 // Periodicities
212 bool m_perx, m_pery;
213 double m_sx, m_sy;
214
215 // Signals
216 int nFourier;
217 std::string m_scellTypeFourier;
218 bool fperx, fpery;
219 int mxmin, mxmax, mymin, mymax;
220 int mfexp;
221
222 std::vector<std::string> m_readout;
223
224 // Wires
225 unsigned int m_nWires;
226 struct wire {
227 double x, y; //< Location.
228 double d; //< Diameter.
229 double v; //< Potential.
230 double e; //< Charge.
231 std::string type; //< Label.
232 double u; //< Length.
233 int ind; //< Readout group.
234 /// Trap radius. Particle is "trapped" if within nTrap * radius of wire.
235 int nTrap;
236 };
237 std::vector<wire> m_w;
238
239 // Stretching weight
240 std::vector<double> weight;
241 // Density
242 std::vector<double> dens;
243 // Mirror charges for force calculations
244 std::vector<double> cnalso;
245
246 // Option for computation of dipole terms
247 bool dipole;
248 // Dipole angle and amplitude
249 std::vector<double> cosph2;
250 std::vector<double> sinph2;
251 std::vector<double> amp2;
252
253 // Parameters for B2 type cells
254 std::vector<double> m_b2sin;
255 // Parameters for C type cells
256 int m_mode;
257 std::complex<double> m_zmult;
258 double m_p1, m_p2, m_c1;
259 // Parameters for D3 type cells
260 // Conformal mapping in polygons
261 std::vector<std::complex<double> > wmap;
262 double m_kappa;
263 // Tables of coefficients
264 std::vector<std::vector<double> > m_cc1;
265 std::vector<std::vector<double> > m_cc2;
266
267 // Reference potential
268 double m_v0;
269 double m_corvta, m_corvtb, m_corvtc;
270
271 // Planes
272 // Existence
273 bool m_ynplan[4];
274 bool m_ynplax, m_ynplay;
275 // Coordinates
276 double m_coplan[4];
277 double m_coplax, m_coplay;
278 // Voltages
279 double m_vtplan[4];
280
281 struct strip {
282 std::string type; //< Label.
283 int ind; //< Readout group.
284 double smin, smax; //< Coordinates.
285 double gap; //< Distance to the opposite electrode.
286 };
287
288 struct pixel {
289 std::string type; //< Label.
290 int ind; //< Readout group.
291 double smin, smax; //< Coordinates in x/y.
292 double zmin, zmax; //< Coordinates in z.
293 double gap; //< Distance to the opposite electrode.
294 };
295
296 struct plane {
297 std::string type; //< Label.
298 int ind; //< Readout group.
299 double ewxcor, ewycor; //< Background weighting fields
300 std::vector<strip> strips1; //< x/y strips.
301 std::vector<strip> strips2; //< z strips.
302 std::vector<pixel> pixels; //< Pixels.
303 };
304
305 std::vector<plane> planes;
306
307 // Tube
308 bool m_tube;
309 int m_mtube, m_ntube;
310 double m_cotube;
311 double m_cotube2;
312 double m_vttube;
313
314 // Capacitance matrix
315 std::vector<std::vector<double> > m_a;
316 // Signal matrix
317 std::vector<std::vector<std::complex<double> > > m_sigmat;
318 // Induced charges on planes
319 std::vector<std::vector<double> > m_qplane;
320
321 // Point charges
322 struct charge3d {
323 double x, y, z; //< Coordinates.
324 double e; //< Charge.
325 };
326 std::vector<charge3d> m_ch3d;
327 unsigned int m_nTermBessel;
328 unsigned int m_nTermPoly;
329
330 // Gravity
331 double down[3];
332
333 void UpdatePeriodicity();
334 void Reset() { CellInit(); }
335
336 void CellInit();
337 bool Prepare();
338 bool CellCheck();
339 bool CellType();
340 bool PrepareStrips();
341
342 bool PrepareSignals();
343 bool SetupWireSignals();
344 bool SetupPlaneSignals();
345
346 // Calculation of charges
347 bool Setup();
348 bool SetupA00();
349 bool SetupB1X();
350 bool SetupB1Y();
351 bool SetupB2X();
352 bool SetupB2Y();
353 bool SetupC10();
354 bool SetupC2X();
355 bool SetupC2Y();
356 bool SetupC30();
357 bool SetupD10();
358 bool SetupD20();
359 bool SetupD30();
360
361 bool IprA00(const int mx, const int my);
362 bool IprB2X(const int my);
363 bool IprB2Y(const int mx);
364 bool IprC2X();
365 bool IprC2Y();
366 bool IprC30();
367 bool IprD10();
368 bool IprD30();
369
370 bool SetupDipole() { return true; }
371
372 // Inversion of capacitance matrix
373 bool Charge();
374
375 // Evaluation of the electric field
376 int Field(const double xin, const double yin, const double zin, double& ex,
377 double& ey, double& ez, double& volt, const bool opt);
378 void FieldA00(const double xpos, const double ypos, double& ex, double& ey,
379 double& volt, const bool opt) const;
380 void FieldB1X(const double xpos, const double ypos, double& ex, double& ey,
381 double& volt, const bool opt) const;
382 void FieldB1Y(const double xpos, const double ypos, double& ex, double& ey,
383 double& volt, const bool opt) const;
384 void FieldB2X(const double xpos, const double ypos, double& ex, double& ey,
385 double& volt, const bool opt) const;
386 void FieldB2Y(const double xpos, const double ypos, double& ex, double& ey,
387 double& volt, const bool opt) const;
388 void FieldC10(const double xpos, const double ypos, double& ex, double& ey,
389 double& volt, const bool opt) const;
390 void FieldC2X(const double xpos, const double ypos, double& ex, double& ey,
391 double& volt, const bool opt) const;
392 void FieldC2Y(const double xpos, const double ypos, double& ex, double& ey,
393 double& volt, const bool opt) const;
394 void FieldC30(const double xpos, const double ypos, double& ex, double& ey,
395 double& volt, const bool opt) const;
396 void FieldD10(const double xpos, const double ypos, double& ex, double& ey,
397 double& volt, const bool opt) const;
398 void FieldD20(const double xpos, const double ypos, double& ex, double& ey,
399 double& volt, const bool opt) const;
400 void FieldD30(const double xpos, const double ypos, double& ex, double& ey,
401 double& volt, const bool opt) const;
402
403 // Field due to point charges
404 void Field3dA00(const double x, const double y, const double z, double& ex,
405 double& ey, double& ez, double& volt);
406 void Field3dB2X(const double x, const double y, const double z, double& ex,
407 double& ey, double& ez, double& volt);
408 void Field3dB2Y(const double x, const double y, const double z, double& ex,
409 double& ey, double& ez, double& volt);
410 void Field3dD10(const double x, const double y, const double z, double& ex,
411 double& ey, double& ez, double& volt);
412 // Evaluation of the weighting field
413 bool Wfield(const double xpos, const double ypos, const double zpos,
414 double& ex, double& ey, double& ez, double& volt,
415 const std::string& label, const bool opt) const;
416 void WfieldWireA00(const double xpos, const double ypos, double& ex,
417 double& ey, double& volt, const int mx, const int my,
418 const int sw, const bool opt) const;
419 void WfieldWireB2X(const double xpos, const double ypos, double& ex,
420 double& ey, double& volt, const int my, const int sw,
421 const bool opt) const;
422 void WfieldWireB2Y(const double xpos, const double ypos, double& ex,
423 double& ey, double& volt, const int mx, const int sw,
424 const bool opt) const;
425 void WfieldWireC2X(const double xpos, const double ypos, double& ex,
426 double& ey, double& volt, const int sw,
427 const bool opt) const;
428 void WfieldWireC2Y(const double xpos, const double ypos, double& ex,
429 double& ey, double& volt, const int sw,
430 const bool opt) const;
431 void WfieldWireC30(const double xpos, const double ypos, double& ex,
432 double& ey, double& volt, const int sw,
433 const bool opt) const;
434 void WfieldWireD10(const double xpos, const double ypos, double& ex,
435 double& ey, double& volt, const int sw,
436 const bool opt) const;
437 void WfieldWireD30(const double xpos, const double ypos, double& ex,
438 double& ey, double& volt, const int sw,
439 const bool opt) const;
440 void WfieldPlaneA00(const double xpos, const double ypos, double& ex,
441 double& ey, double& volt, const int mx, const int my,
442 const int iplane, const bool opt) const;
443 void WfieldPlaneB2X(const double xpos, const double ypos, double& ex,
444 double& ey, double& volt, const int my, const int iplane,
445 const bool opt) const;
446 void WfieldPlaneB2Y(const double xpos, const double ypos, double& ex,
447 double& ey, double& volt, const int mx, const int iplane,
448 const bool opt) const;
449 void WfieldPlaneC2X(const double xpos, const double ypos, double& ex,
450 double& ey, double& volt, const int iplane,
451 const bool opt) const;
452 void WfieldPlaneC2Y(const double xpos, const double ypos, double& ex,
453 double& ey, double& volt, const int iplane,
454 const bool opt) const;
455 void WfieldPlaneC30(const double xpos, const double ypos, double& ex,
456 double& ey, double& volt, const int iplane,
457 const bool opt) const;
458 void WfieldPlaneD10(const double xpos, const double ypos, double& ex,
459 double& ey, double& volt, const int iplane,
460 const bool opt) const;
461 void WfieldPlaneD30(const double xpos, const double ypos, double& ex,
462 double& ey, double& volt, const int iplane,
463 const bool opt) const;
464 void WfieldStripZ(const double xpos, const double ypos, double& ex,
465 double& ey, double& volt, const int ip, const int is,
466 const bool opt) const;
467 void WfieldStripXy(const double xpos, const double ypos, const double zpos,
468 double& ex, double& ey, double& ez, double& volt,
469 const int ip, const int is, const bool opt) const;
470 void WfieldPixel(const double xpos, const double ypos, const double zpos,
471 double& ex, double& ey, double& ez, double& volt,
472 const int ip, const int is, const bool opt) const;
473
474 // Auxiliary functions for C type cells
475 double Ph2(const double xpos, const double ypos) const;
476 double Ph2Lim(const double radius) const {
477 return -log(abs(m_zmult) * radius * (1. - 3. * m_p1 + 5. * m_p2));
478 }
479 void E2Sum(const double xpos, const double ypos,
480 double& ex, double& ey) const;
481
482 // Mapping function for D30 type cells
483 void ConformalMap(const std::complex<double>& z, std::complex<double>& ww,
484 std::complex<double>& wd) const;
485 void InitializeCoefficientTables();
486
487 bool InTube(const double x0, const double y0, const double a,
488 const int n) const;
489
490 // Transformation between cartesian and polar coordinates
491 void Cartesian2Polar(const double x0, const double y0, double& r,
492 double& theta) {
493
494 if (x0 == 0. && y0 == 0.) {
495 r = theta = 0.;
496 return;
497 }
498 r = sqrt(x0 * x0 + y0 * y0);
499 theta = atan2(y0, x0) * RadToDegree;
500 }
501
502 void Polar2Cartesian(const double r, const double theta,
503 double& x0, double& y0) const {
504
505 const double thetap = theta * DegreeToRad;
506 x0 = r * cos(thetap);
507 y0 = r * sin(thetap);
508 }
509
510 // Transformation (r, theta) to (rho, phi) via the map
511 // (r, theta) = (exp(rho), 180 * phi / Pi).
512 void RTheta2RhoPhi(const double rho, const double phi, double& r,
513 double& theta) const {
514
515 r = exp(rho);
516 theta = RadToDegree * phi;
517 }
518};
519}
520
521#endif
bool GetVoltageRange(double &pmin, double &pmax)
Calculate the voltage range [V].
double WeightingPotential(const double x, const double y, const double z, const std::string &label)
void SetPeriodicityX(const double s)
Set the periodic length [cm] in the x-direction.
void ElectricField(const double x, const double y, const double z, double &ex, double &ey, double &ez, double &v, Medium *&m, int &status)
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 EnableChargeCheck(const bool on=true)
bool IsInTrapRadius(const double q0, const double x0, const double y0, const double z0, double &xw, double &yx, double &rw)
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.)
void AddCharge(const double x, const double y, const double z, const double q)
Add a point charge.
void WeightingField(const double x, const double y, const double z, double &wx, double &wy, double &wz, const std::string &label)
bool GetWire(const unsigned int i, double &x, double &y, double &diameter, double &voltage, std::string &label, double &length, double &charge, int &ntrap) const
void AddStripOnPlaneX(const char direction, const double x, const double smin, const double smax, const std::string &label, const double gap=-1.)
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) .
bool GetTube(double &r, double &voltage, int &nEdges, std::string &label) const
bool GetPlaneX(const unsigned int i, double &x, double &voltage, std::string &label) const
void AddPlaneY(const double y, const double voltage, const std::string &label)
Add a plane at constant y.
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.)
bool GetBoundingBox(double &x0, double &y0, double &z0, double &x1, double &y1, double &z1)
Get the bounding box coordinates.
void AddStripOnPlaneY(const char direction, const double y, const double smin, const double smax, const std::string &label, const double gap=-1.)
void AddReadout(const std::string &label)
Setup the weighting field for a given group of wires or planes.
void ElectricField(const double x, const double y, const double z, double &ex, double &ey, double &ez, Medium *&m, int &status)
void SetPeriodicityY(const double s)
Set the periodic length [cm] in the y-direction.
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)
bool GetPlaneY(const unsigned int i, double &y, double &voltage, std::string &label) const
Abstract base class for components.
virtual Medium * GetMedium(const double x, const double y, const double z)
Get the medium at a given location (x, y, z).
Abstract base class for media.
Definition: Medium.hh:11
bool IsDriftable() const
Definition: Medium.hh:57
DoubleAc cos(const DoubleAc &f)
Definition: DoubleAc.cpp:432
DoubleAc exp(const DoubleAc &f)
Definition: DoubleAc.cpp:377
DoubleAc sin(const DoubleAc &f)
Definition: DoubleAc.cpp:384
DoubleAc sqrt(const DoubleAc &f)
Definition: DoubleAc.cpp:314