Garfield++ 4.0
A toolkit for the detailed simulation of particle detectors based on ionisation measurement in gases and semiconductors
Loading...
Searching...
No Matches
ComponentParallelPlate.cc
Go to the documentation of this file.
1#include "Garfield/ComponentParallelPlate.hh"
2
3#include <TF1.h>
4#include <TF2.h>
5
6#include <cmath>
7#include <limits>
8#include <iostream>
9
10#include "Garfield/GarfieldConstants.hh"
11
12namespace Garfield {
13
14ComponentParallelPlate::ComponentParallelPlate() : Component("ParallelPlate") {}
15
16void ComponentParallelPlate::Setup(double g, double b, double eps, double V,
17 double sigma) {
18 // TODO: can g, b be negative?
19 m_g = g;
20 m_b = b;
21 if (eps < 1.) {
22 std::cerr << m_className << "::Setup: Epsilon must be >= 1.\n";
23 return;
24 }
25 m_eps = eps;
26 m_V = V;
27 // TODO: can sigma be negative?
28 m_sigma = sigma;
29 if (sigma == 0) {
30 m_ezg = -m_eps * m_V / (m_b + m_eps * m_g);
31 m_ezb = -m_V / (m_b + m_eps * m_g);
32 } else {
33 // For large times the resistive layer will act as a perfect conductor.
34 m_ezg = -m_V / m_g;
35 m_ezb = 0.;
36 }
37 std::cout << m_className << "::Setup: Geometry set.\n";
38}
39
40bool ComponentParallelPlate::GetBoundingBox(double& x0, double& y0, double& z0,
41 double& x1, double& y1,
42 double& z1) {
43 // If a geometry is present, try to get the bounding box from there.
44 if (m_geometry) {
45 if (m_geometry->GetBoundingBox(x0, y0, z0, x1, y1, z1)) return true;
46 }
47 x0 = -std::numeric_limits<double>::infinity();
48 y0 = -std::numeric_limits<double>::infinity();
49 x1 = +std::numeric_limits<double>::infinity();
50 y1 = +std::numeric_limits<double>::infinity();
51 // TODO: check!
52 z0 = 0.;
53 z1 = m_g;
54 return true;
55}
56
57double ComponentParallelPlate::IntegrateField(const Electrode& el, int comp,
58 const double x, const double y,
59 const double z) {
60 switch (el.ind) {
61 case structureelectrode::Plane: {
62 if (comp == fieldcomponent::zcomp)
63 return m_eps * m_Vw / (m_b + m_eps * m_g);
64 return 0.;
65 break;
66 }
67 case structureelectrode::Pixel: {
68 auto WFieldPixel = [=](double* k, double* /*p*/) {
69 double kx = k[0];
70 double ky = k[1];
71
72 double K = std::sqrt(kx * kx + ky * ky);
73
74 double intsol = 1.;
75
76 switch (comp) {
77 case fieldcomponent::xcomp: {
78 intsol *= 1. / (ky * cosh(m_g * K) * sinh(m_b * K) +
79 m_eps * ky * cosh(m_b * K) * sinh(m_g * K));
80
81 intsol *= cos(ky * (y - el.ypos)) * sin((kx * el.lx) / 2) *
82 sin((ky * el.ly) / 2) * sin(kx * (x - el.xpos)) *
83 sinh(K * (m_g - z));
84 break;
85 }
86 case fieldcomponent::ycomp: {
87 intsol *= 1. / (kx * cosh(m_g * K) * sinh(m_b * K) +
88 m_eps * kx * cosh(m_b * K) * sinh(m_g * K));
89 intsol *= sin(ky * (y - el.ypos)) * sin((kx * el.lx) / 2) *
90 sin((ky * el.ly) / 2) * cos(kx * (x - el.xpos)) *
91 sinh(K * (m_g - z));
92 break;
93 }
94 case fieldcomponent::zcomp: {
95 intsol *= 1. / (ky * kx * cosh(m_g * K) * sinh(m_b * K) +
96 m_eps * ky * kx * cosh(m_b * K) * sinh(m_g * K));
97 intsol *= K * cos(ky * (y - el.ypos)) * sin((kx * el.lx) / 2) *
98 sin((ky * el.ly) / 2) * cos(kx * (x - el.xpos)) *
99 cosh(K * (m_g - z));
100 break;
101 }
102 }
103 return intsol;
104 };
105 TF2* fw =
106 new TF2("WFieldPixel", WFieldPixel, 0, 10 * m_g, 0, 10 * m_g, 0);
107 const double sol = fw->Integral(0, 10 * m_g, 0, 10 * m_g, 1.e-6);
108 delete fw;
109 return (4 * m_eps * m_Vw / Pi2) * sol;
110
111 break;
112 }
113 case structureelectrode::Strip: {
114 auto WFieldStrip = [=](double* k, double* /*p*/) {
115 double kk = k[0];
116
117 double intsol = 1. / ((cosh(m_g * kk) * sinh(m_b * kk) +
118 m_eps * cosh(m_b * kk) * sinh(m_g * kk)));
119 switch (comp) {
120 case fieldcomponent::xcomp: {
121 intsol *= (sin(kk * el.lx / 2) * sin(kk * (x - el.xpos)) *
122 sinh(kk * (m_g - z)));
123 break;
124 }
125 case fieldcomponent::zcomp: {
126 intsol *= (sin(kk * el.lx / 2) * cos(kk * (x - el.xpos)) *
127 cosh(kk * (m_g - z)));
128 break;
129 }
130 }
131
132 return intsol;
133 };
134
135 if (comp == ycomp) {
136 return 0.;
137 }
138 TF1* fw = new TF1("WFieldStrip", WFieldStrip, 0, 10 * m_g, 0);
139 double sol = fw->Integral(0, 10 * m_g);
140 delete fw;
141 return (2 * m_eps * m_Vw / Pi) * sol;
142 break;
143 }
144 default: {
145 std::cerr << m_className << "::IntegrateField: Unknown electrode type.\n";
146 return 0.;
147 }
148 }
149}
150
151double ComponentParallelPlate::IntegrateDelayedField(const Electrode& el,
152 int comp, const double x,
153 const double y,
154 const double z,
155 const double t) {
156 switch (el.ind) {
157 case structureelectrode::Plane: {
158 if (comp == fieldcomponent::zcomp)
159 return m_eps * m_Vw *
160 (1 - exp(-t * m_g * m_sigma / (m_eps0 * (m_b + m_eps * m_g)))) /
161 (m_b + m_eps * m_g);
162 return 0.;
163 break;
164 }
165 case structureelectrode::Pixel: {
166 auto WFieldPixel = [=](double* k, double* /*p*/) {
167 double kx = k[0];
168 double ky = k[1];
169
170 double K = std::sqrt(kx * kx + ky * ky);
171
172 double tau = m_eps0 *
173 (m_eps + cosh(m_g * K) * sinh(m_b * K) /
174 (cosh(m_b * K) * sinh(m_g * K))) *
175 (1 / m_sigma);
176
177 double intsol = 1. / (cosh(m_g * K) * sinh(m_b * K) +
178 m_eps * cosh(m_b * K) * sinh(m_g * K));
179
180 switch (comp) {
181 case fieldcomponent::xcomp: {
182 intsol *= (1 - exp(-t / tau)) * cos(ky * (y - el.ypos)) *
183 cosh(m_g * K) * sin((kx * el.lx) / 2) *
184 sin(kx * (x - el.xpos)) * sinh(K * (m_g - z)) *
185 tanh(m_b * K) / (ky * sinh(m_g * K));
186 break;
187 }
188 case fieldcomponent::ycomp: {
189 intsol *= (1 - exp(-t / tau)) * sin(ky * (y - el.ypos)) *
190 cosh(m_g * K) * sin((kx * el.lx) / 2) *
191 cos(kx * (x - el.xpos)) * cosh(K * (m_g - z)) *
192 tanh(m_b * K) / (kx * sinh(m_g * K));
193 break;
194 }
195 case fieldcomponent::zcomp: {
196 intsol *= (1 - exp(-t / tau)) * cos(ky * (y - el.ypos)) *
197 cosh(m_g * K) * sin((K * el.lx) / 2) *
198 cos(kx * (x - el.xpos)) * cosh(K * (m_g - z)) *
199 tanh(m_b * K) / (kx * ky * sinh(m_g * K));
200 break;
201 }
202 }
203 return intsol;
204 };
205
206 TF2* fw =
207 new TF2("WFieldPixel", WFieldPixel, 0, 10 * m_g, 0, 10 * m_g, 0);
208 const double sol = fw->Integral(0, 10 * m_g, 0, 10 * m_g, 1.e-6);
209 delete fw;
210 return (4 * m_eps * m_Vw / Pi2) * sol;
211 break;
212 }
213 case structureelectrode::Strip: {
214 auto WFieldStrip = [=](double* k, double* /*p*/) {
215 double kk = k[0];
216 double tau = m_eps0 *
217 (m_eps + cosh(m_g * kk) * sinh(m_b * kk) /
218 (cosh(m_b * kk) * sinh(m_g * kk))) *
219 (1 / m_sigma);
220
221 double intsol = 1 / (cosh(m_g * kk) * sinh(m_b * kk) +
222 m_eps * cosh(m_b * kk) * sinh(m_g * kk));
223 switch (comp) {
224 case fieldcomponent::xcomp: {
225 intsol *= (1 - exp(-t / tau)) * cosh(m_g * kk) *
226 sin((kk * el.lx) / 2) * sin(kk * (x - el.xpos)) *
227 sinh(kk * (m_g - z)) * tanh(m_b * kk) / sinh(m_g * kk);
228 break;
229 }
230 case fieldcomponent::zcomp: {
231 intsol *= (1 - exp(-t / tau)) * cosh(m_g * kk) *
232 sin((kk * el.lx) / 2) * cos(kk * (x - el.xpos)) *
233 cosh(kk * (m_g - z)) * tanh(m_b * kk) / sinh(m_g * kk);
234 break;
235 }
236 }
237
238 return intsol;
239 };
240
241 if (comp == ycomp) {
242 return 0.;
243 }
244 TF1* fw = new TF1("WFieldStrip", WFieldStrip, 0, 10 * m_g, 0);
245 const double sol = fw->Integral(0, 10 * m_g);
246 delete fw;
247 return (2 * m_eps * m_Vw / Pi) * sol;
248 break;
249 }
250 default: {
251 std::cerr << m_className << "::IntegrateDelayedField:\n"
252 << " Unknown electrode type.\n";
253 return 0.;
254 }
255 }
256}
257
258double ComponentParallelPlate::IntegratePromptPotential(const Electrode& el,
259 const double x,
260 const double y,
261 const double z) {
262 switch (el.ind) {
263 case structureelectrode::Plane: {
264 double sol = m_eps * m_Vw * (m_g - z) / (m_b + m_eps * m_g);
265 return std::abs(sol) > m_precision ? sol : 0.;
266 break;
267 }
268 case structureelectrode::Pixel: {
269 auto WPotentialPixel = [=](double* k, double* /*p*/) {
270 double kx = k[0];
271 double ky = k[1];
272
273 double K = std::sqrt(kx * kx + ky * ky);
274
275 double intsol = 1.;
276
277 intsol *= cos(kx * (x - el.xpos)) * sin(kx * el.lx / 2) *
278 cos(ky * (y - el.ypos)) * sin(ky * el.ly / 2) *
279 sinh(K * (m_g - z)) /
280 (kx * ky *
281 (sinh(m_b * K) * cosh(m_g * K) +
282 m_eps * sinh(m_g * K) * cosh(m_b * K)));
283
284 return intsol;
285 };
286
287 TF2* pw = new TF2("WPotentialPixel", WPotentialPixel, 0, 10 * m_g, 0,
288 10 * m_g, 0);
289 const double sol = pw->Integral(0, 2 * m_g, 0, 2 * m_g, 1.e-6);
290 delete pw;
291 return (4 * m_eps * m_Vw / Pi2) * sol;
292 break;
293 }
294 case structureelectrode::Strip: {
295 auto WPotentialStrip = [=](double* k, double* /*p*/) {
296 double kk = k[0];
297
298 double intsol = 1. / (kk * (cosh(m_g * kk) * sinh(m_b * kk) +
299 m_eps * cosh(m_b * kk) * sinh(m_g * kk)));
300 intsol *= (sin(kk * el.lx / 2) * cos(kk * (x - el.xpos)) *
301 sinh(kk * (m_g - z)));
302
303 return intsol;
304 };
305
306 TF1* pw = new TF1("WPotentialStrip", WPotentialStrip, 0, 10 * m_g, 0);
307 const double sol = pw->Integral(0, 10 * m_g);
308 delete pw;
309 return (2 * m_eps * m_Vw / Pi) * sol;
310 break;
311 }
312 default: {
313 std::cerr << m_className << "::IntegratePromptPotential:\n"
314 << " Unknown electrode type.\n";
315 return 0.;
316 }
317 }
318}
319
320double ComponentParallelPlate::IntegrateDelayedPotential(const Electrode& el,
321 const double x,
322 const double y,
323 const double z,
324 const double t) {
325 switch (el.ind) {
326 case structureelectrode::Plane: {
327 double tau =
328 m_eps0 * (m_eps + m_b / m_g) /
329 (m_sigma); // Note to self: You dropt the eps) here for convenience.
330
331 double sol = m_Vw * (1 - exp(-t / tau)) * (m_b * (m_g - z)) /
332 (m_g * (m_b + m_eps * m_g));
333 return std::abs(sol) > m_precision ? sol : 0.;
334 break;
335 }
336 case structureelectrode::Pixel: {
337 auto WPotentialPixel = [=](double* k, double* /*p*/) {
338 double kx = k[0];
339 double ky = k[1];
340
341 double K = std::sqrt(kx * kx + ky * ky);
342 double tau = m_eps0 *
343 (m_eps + cosh(m_g * K) * sinh(m_b * K) /
344 (cosh(m_b * K) * sinh(m_g * K))) *
345 (1 / m_sigma); // Note to self: You dropt the eps) here
346 // for convenience.
347
348 double intsol = 1. / (kx * ky *
349 (sinh(m_b * K) * cosh(m_g * K) +
350 m_eps * sinh(m_g * K) * cosh(m_b * K)));
351
352 intsol *= cos(kx * (x - el.xpos)) * sin(kx * el.lx / 2) *
353 cos(ky * (y - el.ypos)) * sin(ky * el.ly / 2) *
354 sinh(K * (m_g - z)) * tanh(m_b * K) * cosh(m_g * K) *
355 (1 - exp(-t / tau)) / sinh(m_g * K);
356
357 return intsol;
358 };
359
360 TF2* pw = new TF2("WPotentialPixel", WPotentialPixel, 0, 10 * m_g, 0,
361 10 * m_g, 0);
362 const double sol = pw->Integral(0, 2 * m_g, 0, 2 * m_g, 1.e-20);
363 delete pw;
364 return (4 * m_Vw / Pi2) * sol;
365 break;
366 }
367 case structureelectrode::Strip: {
368 auto WPotentialStrip = [=](double* k, double* /*p*/) {
369 double kk = k[0];
370
371 double tau = m_eps0 *
372 (m_eps + cosh(m_g * kk) * sinh(m_b * kk) /
373 (cosh(m_b * kk) * sinh(m_g * kk))) *
374 (1 / m_sigma);
375
376 double intsol = 1. / (kk * (cosh(m_g * kk) * sinh(m_b * kk) +
377 m_eps * cosh(m_b * kk) * sinh(m_g * kk)));
378 intsol *= (sin(kk * el.lx / 2) * cos(kk * (x - el.xpos)) *
379 sinh(kk * (m_g - z)) * cosh(m_g * kk) * tanh(m_b * kk)) *
380 (1 - exp(-t / tau)) / sinh(m_g * kk);
381
382 return intsol;
383 };
384
385 TF1* pw = new TF1("WPotentialStrip", WPotentialStrip, 0, 10 * m_g, 0);
386 const double sol = pw->Integral(0, 8 * m_g);
387 delete pw;
388 return (2 * m_Vw / Pi) * sol;
389 break;
390 }
391 default: {
392 std::cerr << m_className << "::IntegrateDelayedPotential:\n"
393 << " Unknown electrode type.\n";
394 return 0.;
395 }
396 }
397}
398
399void ComponentParallelPlate::ElectricField(const double x, const double y,
400 const double z, double& ex,
401 double& ey, double& ez, Medium*& m,
402 int& status) {
403 ex = ey = 0.;
404
405 if (z < 0) {
406 ez = m_ezb;
407 } else {
408 ez = m_ezb;
409 }
410
411 m = m_geometry ? m_geometry->GetMedium(x, y, z) : m_medium;
412
413 if (!m) {
414 if (m_debug) {
415 std::cout << m_className << "::ElectricField: No medium at (" << x << ", "
416 << y << ", " << z << ").\n";
417 }
418 status = -6;
419 return;
420 }
421
422 if (z > 0) {
423 status = 0;
424 } else {
425 status = -5;
426 }
427}
428
429void ComponentParallelPlate::ElectricField(const double x, const double y,
430 const double z, double& ex,
431 double& ey, double& ez, double& v,
432 Medium*& m, int& status) {
433 ex = ey = 0.;
434
435 if (z > 0.) {
436 ez = m_ezg;
437 } else {
438 ez = m_ezb;
439 }
440
441 if (m_sigma == 0) {
442 v = -m_eps * m_V * (m_g - z) / (m_b + m_eps * m_g);
443 } else {
444 v = -m_eps * m_V * (m_g - z) / (m_eps * m_g);
445 }
446
447 m = m_geometry ? m_geometry->GetMedium(x, y, z) : m_medium;
448 if (!m) {
449 if (m_debug) {
450 std::cout << m_className << "::ElectricField: No medium at (" << x << ", "
451 << y << ", " << z << ").\n";
452 }
453 status = -6;
454 return;
455 }
456
457 if (z > 0) {
458 status = 0;
459 } else {
460 status = -5;
461 }
462}
463
464bool ComponentParallelPlate::GetVoltageRange(double& vmin, double& vmax) {
465 if (m_V == 0) return false;
466
467 if (m_V < 0) {
468 vmin = m_V;
469 vmax = 0;
470 } else {
471 vmin = 0;
472 vmax = m_V;
473 }
474 return true;
475}
476
477void ComponentParallelPlate::WeightingField(const double x, const double y,
478 const double z, double& wx,
479 double& wy, double& wz,
480 const std::string& label) {
481 wx = 0;
482 wy = 0;
483 wz = 0;
484
485 for (const auto& electrode : m_readout_p) {
486 if (electrode.label == label) {
487 wx = electrode.flip *
488 IntegrateField(electrode, fieldcomponent::xcomp, x, y, z);
489 wy = electrode.flip *
490 IntegrateField(electrode, fieldcomponent::ycomp, x, y, z);
491 wz = electrode.flip *
492 IntegrateField(electrode, fieldcomponent::zcomp, x, y, z);
493 }
494 }
495}
496
497double ComponentParallelPlate::WeightingPotential(const double x,
498 const double y,
499 const double z,
500 const std::string& label) {
501 double ret = 0.;
502
503 for (const auto& electrode : m_readout_p) {
504 if (electrode.label == label) {
505 if (!electrode.m_usegrid) {
506 ret += electrode.flip * IntegratePromptPotential(electrode, x, y, z);
507 } else {
508 ret += FindWeightingPotentialInGrid(electrode, x, y, z);
509 }
510 }
511 }
512 return ret;
513}
514
515double ComponentParallelPlate::DelayedWeightingPotential(
516 const double x, const double y, const double z, const double t,
517 const std::string& label) {
518 if (m_sigma == 0) {
519 if (m_debug) {
520 std::cout << m_className << "::DelayedWeightingPotential:\n"
521 << " Conductivity is set to zero.\n";
522 }
523 return 0.;
524 }
525
526 double ret = 0.;
527
528 for (const auto& electrode : m_readout_p) {
529 if (electrode.label == label) {
530 if (!electrode.m_usegrid) {
531 ret +=
532 electrode.flip * IntegrateDelayedPotential(electrode, x, y, z, t);
533 } else {
534 ret += FindDelayedWeightingPotentialInGrid(electrode, x, y, z, t);
535 }
536 }
537 }
538
539 return ret;
540}
541
542void ComponentParallelPlate::DelayedWeightingField(
543 const double x, const double y, const double z, const double t, double& wx,
544 double& wy, double& wz, const std::string& label) {
545 wx = 0.;
546 wy = 0.;
547 wz = 0.;
548
549 if (m_sigma == 0) {
550 if (m_debug) {
551 std::cout << m_className << "::DelayedWeightingField:\n"
552 << " Conductivity is set to zero.\n";
553 }
554 return;
555 }
556
557 for (const auto& electrode : m_readout_p) {
558 if (electrode.label == label) {
559 wx = electrode.flip *
560 IntegrateDelayedField(electrode, fieldcomponent::xcomp, x, y, z, t);
561 wy = electrode.flip *
562 IntegrateDelayedField(electrode, fieldcomponent::ycomp, x, y, z, t);
563 wz = electrode.flip *
564 IntegrateDelayedField(electrode, fieldcomponent::zcomp, x, y, z, t);
565 }
566 }
567}
568
569void ComponentParallelPlate::Reset() {
570 m_readout.clear();
571 m_readout_p.clear();
572
573 m_g = 0.;
574 m_b = 0.;
575 m_eps = 1.;
576 m_V = 0.;
577}
578
579void ComponentParallelPlate::UpdatePeriodicity() {
580 if (m_debug) {
581 std::cerr << m_className << "::UpdatePeriodicity:\n"
582 << " Periodicities are not supported.\n";
583 }
584}
585
586void ComponentParallelPlate::AddPixel(double x, double y, double lx_input,
587 double ly_input,
588 const std::string& label) {
589 const auto it = std::find(m_readout.cbegin(), m_readout.cend(), label);
590 if (it == m_readout.end() && m_readout.size() > 0) {
591 std::cerr << m_className << "::AddPixel:\n"
592 << "Note that the label " << label << " is already in use.\n";
593 }
594 Electrode pixel;
595 pixel.label = label;
596 pixel.ind = structureelectrode::Pixel;
597 pixel.xpos = x;
598 pixel.ypos = y;
599 pixel.lx = lx_input;
600 pixel.ly = ly_input;
601
602 m_readout.push_back(label);
603 m_readout_p.push_back(std::move(pixel));
604 std::cout << m_className << "::AddPixel: Added pixel electrode.\n";
605}
606
607void ComponentParallelPlate::AddStrip(double x, double lx_input,
608 const std::string& label) {
609 const auto it = std::find(m_readout.cbegin(), m_readout.cend(), label);
610 if (it == m_readout.end() && m_readout.size() > 0) {
611 std::cerr << m_className << "::AddStrip:\n"
612 << "Note that the label " << label << " is already in use.\n";
613 }
614 Electrode strip;
615 strip.label = label;
616 strip.ind = structureelectrode::Strip;
617 strip.xpos = x;
618 strip.lx = lx_input;
619
620 m_readout.push_back(label);
621 m_readout_p.push_back(std::move(strip));
622
623 std::cout << m_className << "::AddStrip: Added strip electrode.\n";
624}
625
626void ComponentParallelPlate::AddPlane(const std::string& label, bool anode) {
627 const auto it = std::find(m_readout.cbegin(), m_readout.cend(), label);
628 if (it == m_readout.end() && m_readout.size() > 0) {
629 std::cerr << m_className << "::AddPlane:\n"
630 << "Note that the label " << label << " is already in use.\n";
631 }
632 Electrode plate;
633 plate.label = label;
634 plate.ind = structureelectrode::Plane;
635
636 if (!anode) plate.flip = -1;
637
638 m_readout.push_back(label);
639 m_readout_p.push_back(std::move(plate));
640
641 std::cout << m_className << "::AddPlane: Added plane electrode.\n";
642}
643
644Medium* ComponentParallelPlate::GetMedium(const double x, const double y,
645 const double z) {
646 if (m_geometry) {
647 return m_geometry->GetMedium(x, y, z);
648 } else if (m_medium) {
649 return m_medium;
650 }
651 return nullptr;
652}
653
654void ComponentParallelPlate::SetWeightingPotentialGrid(
655 const std::string& label, const double xmin, const double xmax,
656 const double xsteps, const double ymin, const double ymax,
657 const double ysteps, const double zmin, const double zmax,
658 const double zsteps, const double tmin, const double tmax,
659 const double tsteps) {
660 for (auto& electrode : m_readout_p) {
661 if (electrode.label == label) {
662 electrode.gridXSteps = xsteps;
663 electrode.gridYSteps = ysteps;
664 electrode.gridZSteps = zsteps;
665 electrode.gridTSteps = tsteps;
666
667 if (xsteps == 0) electrode.gridXSteps = 1;
668 if (ysteps == 0) electrode.gridYSteps = 1;
669
670 electrode.gridX0 = xmin;
671 electrode.gridY0 = ymin;
672 electrode.gridZ0 = zmin;
673 electrode.gridT0 = tmin;
674
675 electrode.gridXStepSize = (xmax - xmin) / xsteps;
676 electrode.gridYStepSize = (ymax - ymin) / ysteps;
677 electrode.gridZStepSize = (zmax - zmin) / zsteps;
678 electrode.gridTStepSize = (tmax - tmin) / tsteps;
679
680 std::vector<double> nhz(zsteps, 0);
681 std::vector<std::vector<double>> nhy(ysteps, nhz);
682 std::vector<std::vector<std::vector<double>>> nhx(xsteps, nhy);
683 electrode.gridPromptV = nhx;
684
685 std::vector<double> nht(tsteps, 0);
686 std::vector<std::vector<double>> nhzd(zsteps, nht);
687 std::vector<std::vector<std::vector<double>>> nhyd(ysteps, nhzd);
688 std::vector<std::vector<std::vector<std::vector<double>>>> nhxd(xsteps,
689 nhyd);
690 electrode.gridDelayedV = nhxd;
691
692 for (int ix = 0; ix < xsteps; ix++) {
693 for (int iy = 0; iy < xsteps; iy++) {
694 for (int iz = 0; iz < xsteps; iz++) {
695 if (iz * zsteps + zmin >= 0)
696 electrode.gridPromptV[ix][iy][iz] =
697 electrode.flip * IntegratePromptPotential(
698 electrode, ix * xsteps + xmin,
699 iy * ysteps + ymin, iz * zsteps + zmin);
700
701 for (int it = 0; it < tsteps; it++) {
702 if (iz * zsteps + zmin >= 0)
703 electrode.gridDelayedV[ix][iy][iz][it] =
704 electrode.flip * IntegrateDelayedPotential(
705 electrode, ix * xsteps + xmin,
706 iy * ysteps + ymin, iz * zsteps + zmin,
707 it * tsteps + tmin);
708 }
709 }
710 }
711 }
712
713 electrode.m_usegrid = true;
714 }
715 }
716}
717
718void ComponentParallelPlate::SetWeightingPotentialGrids(
719 const double xmin, const double xmax, const double xsteps,
720 const double ymin, const double ymax, const double ysteps,
721 const double zmin, const double zmax, const double zsteps,
722 const double tmin, const double tmax, const double tsteps) {
723 for (const auto& electrode : m_readout_p) {
724 SetWeightingPotentialGrid(electrode.label, xmin, xmax, xsteps, ymin, ymax,
725 ysteps, zmin, zmax, zsteps, tmin, tmax, tsteps);
726 }
727}
728
729double ComponentParallelPlate::FindWeightingPotentialInGrid(const Electrode& el,
730 const double x,
731 const double y,
732 const double z) {
733 switch (el.ind) {
734 case structureelectrode::Plane: {
735 return el.flip * IntegratePromptPotential(el, x, y, z);
736 break;
737 }
738 case structureelectrode::Strip: {
739 int ix = floor((x - el.gridX0) / el.gridXStepSize);
740 int iz = floor((z - el.gridZ0) / el.gridZStepSize);
741
742 if (ix < 0 || ix >= el.gridXSteps || iz < 0 || iz >= el.gridZSteps)
743 return IntegratePromptPotential(el, x, y, z);
744
745 double ret = 0;
746
747 for (int i = 0; i < 2; i++) {
748 for (int j = 0; j < 2; j++) {
749 ret += FindWeightFactor(
750 el, std::abs((ix + i) * el.gridXStepSize + el.gridX0 - x), 0,
751 std::abs((iz + j) * el.gridZStepSize + el.gridZ0 - z)) *
752 el.gridPromptV[ix + i][0][iz + j];
753 }
754 }
755
756 return ret;
757 break;
758 }
759 case structureelectrode::Pixel: {
760 int ix = floor((x - el.gridX0) / el.gridXStepSize);
761 int iy = floor((y - el.gridY0) / el.gridYStepSize);
762 int iz = floor((z - el.gridZ0) / el.gridZStepSize);
763
764 if (ix < 0 || ix >= el.gridXSteps || iz < 0 || iz >= el.gridYSteps ||
765 iz < 0 || iz >= el.gridZSteps)
766 return IntegratePromptPotential(el, x, y, z);
767
768 double ret = 0;
769
770 for (int i = 0; i < 2; i++) {
771 for (int j = 0; j < 2; j++) {
772 for (int k = 0; k < 2; k++) {
773 ret += FindWeightFactor(
774 el, std::abs((ix + i) * el.gridXStepSize + el.gridX0 - x),
775 std::abs((iy + k) * el.gridYStepSize + el.gridY0 - y),
776 std::abs((iz + j) * el.gridZStepSize + el.gridZ0 - z)) *
777 el.gridPromptV[ix + i][iy + k][iz + j];
778 }
779 }
780 }
781 return ret;
782 break;
783 }
784 }
785 return 0.;
786}
787
788double ComponentParallelPlate::FindDelayedWeightingPotentialInGrid(
789 const Electrode& el, const double x, const double y, const double z,
790 const double t) {
791 switch (el.ind) {
792 case structureelectrode::Plane: {
793 return el.flip * IntegrateDelayedPotential(el, x, y, z, t);
794 break;
795 }
796 case structureelectrode::Strip: {
797 int ix = floor((x - el.gridX0) / el.gridXStepSize);
798 int iz = floor((z - el.gridZ0) / el.gridZStepSize);
799 int it = floor((t - el.gridT0) / el.gridTStepSize);
800
801 if (ix < 0 || ix >= el.gridXSteps || iz < 0 || iz >= el.gridZSteps ||
802 it < 0 || it >= el.gridTSteps)
803 return IntegrateDelayedPotential(el, x, y, z, t);
804
805 double ret = 0;
806
807 for (int i = 0; i < 2; i++) {
808 for (int j = 0; j < 2; j++) {
809 for (int l = 0; l < 2; l++) {
810 ret += FindWeightFactor(
811 el, std::abs((ix + i) * el.gridXStepSize + el.gridX0 - x), 0,
812 std::abs((iz + j) * el.gridZStepSize + el.gridZ0 - z),
813 std::abs((it + l) * el.gridTStepSize + el.gridT0 - t)) *
814 el.gridDelayedV[ix + i][0][iz + j][it + l];
815 }
816 }
817 }
818
819 return ret;
820 break;
821 }
822 case structureelectrode::Pixel: {
823 int ix = floor((x - el.gridX0) / el.gridXStepSize);
824 int iy = floor((y - el.gridY0) / el.gridYStepSize);
825 int iz = floor((z - el.gridZ0) / el.gridZStepSize);
826 int it = floor((t - el.gridT0) / el.gridTStepSize);
827
828 if (ix < 0 || ix >= el.gridXSteps || iz < 0 || iz >= el.gridYSteps ||
829 iz < 0 || iz >= el.gridZSteps || it < 0 || it >= el.gridTSteps)
830 return IntegrateDelayedPotential(el, x, y, z, t);
831
832 double ret = 0;
833
834 for (int i = 0; i < 2; i++) {
835 for (int j = 0; j < 2; j++) {
836 for (int k = 0; k < 2; k++) {
837 for (int l = 0; l < 2; l++) {
838 ret += FindWeightFactor(
839 el, std::abs((ix + i) * el.gridXStepSize + el.gridX0 - x),
840 std::abs((iy + k) * el.gridYStepSize + el.gridY0 - y),
841 std::abs((iz + j) * el.gridZStepSize + el.gridZ0 - z),
842 std::abs((it + l) * el.gridTStepSize + el.gridT0 - t)) *
843 el.gridDelayedV[ix + i][iy + k][iz + j][it + l];
844 }
845 }
846 }
847 }
848 return ret;
849 break;
850 }
851 }
852 return 0.;
853}
854
855double ComponentParallelPlate::FindWeightFactor(const Electrode& el,
856 const double dx,
857 const double dy,
858 const double dz) {
859 double fact = 0;
860
861 switch (el.ind) {
862 case structureelectrode::Strip: {
863 fact = (el.gridXStepSize - dx) * (el.gridZStepSize - dz) /
864 (el.gridXStepSize * el.gridZStepSize);
865 break;
866 }
867 case structureelectrode::Pixel: {
868 fact = (el.gridXStepSize - dx) * (el.gridYStepSize - dy) *
869 (el.gridZStepSize - dz) /
870 (el.gridXStepSize * el.gridYStepSize * el.gridZStepSize);
871 break;
872 }
873 }
874
875 return fact;
876}
877
878double ComponentParallelPlate::FindWeightFactor(const Electrode& el,
879 const double dx,
880 const double dy,
881 const double dz,
882 const double dt) {
883 double fact = 0;
884
885 switch (el.ind) {
886 case structureelectrode::Strip: {
887 fact = (el.gridXStepSize - dx) * (el.gridZStepSize - dz) *
888 (el.gridXStepSize - dt) /
889 (el.gridXStepSize * el.gridZStepSize * el.gridTStepSize);
890 break;
891 }
892 case structureelectrode::Pixel: {
893 fact = (el.gridXStepSize - dx) * (el.gridYStepSize - dy) *
894 (el.gridZStepSize - dz) * (el.gridXStepSize - dt) /
895 (el.gridXStepSize * el.gridYStepSize * el.gridZStepSize *
896 el.gridTStepSize);
897 break;
898 }
899 }
900
901 return fact;
902}
903
904} // namespace Garfield
Garfield::ComponentParallelPlate::Setup
void Setup(double g, double b, double eps, double v, double sigma=0.)
Definition: ComponentParallelPlate.cc:16
Garfield::ComponentParallelPlate::DelayedWeightingField
void DelayedWeightingField(const double x, const double y, const double z, const double t, double &wx, double &wy, double &wz, const std::string &label) override
Definition: ComponentParallelPlate.cc:542
Garfield::ComponentParallelPlate::WeightingPotential
double WeightingPotential(const double x, const double y, const double z, const std::string &label) override
Definition: ComponentParallelPlate.cc:497
Garfield::ComponentParallelPlate::AddStrip
void AddStrip(double x, double lx, const std::string &label)
Add strip electrode.
Definition: ComponentParallelPlate.cc:607
Garfield::ComponentParallelPlate::DelayedWeightingPotential
double DelayedWeightingPotential(const double x, const double y, const double z, const double t, const std::string &label) override
Definition: ComponentParallelPlate.cc:515
Garfield::ComponentParallelPlate::WeightingField
void WeightingField(const double x, const double y, const double z, double &wx, double &wy, double &wz, const std::string &label) override
Definition: ComponentParallelPlate.cc:477
Garfield::ComponentParallelPlate::ElectricField
void ElectricField(const double x, const double y, const double z, double &ex, double &ey, double &ez, Medium *&m, int &status) override
Definition: ComponentParallelPlate.cc:399
Garfield::ComponentParallelPlate::AddPlane
void AddPlane(const std::string &label, bool anode=true)
Definition: ComponentParallelPlate.cc:626
Garfield::ComponentParallelPlate::SetWeightingPotentialGrid
void SetWeightingPotentialGrid(const std::string &label, const double xmin, const double xmax, const double xsteps, const double ymin, const double ymax, const double ysteps, const double zmin, const double zmax, const double zsteps, const double tmin, const double tmax, const double tsteps)
Definition: ComponentParallelPlate.cc:654
Garfield::ComponentParallelPlate::GetMedium
Medium * GetMedium(const double x, const double y, const double z) override
Get the medium at a given location (x, y, z).
Definition: ComponentParallelPlate.cc:644
Garfield::ComponentParallelPlate::SetWeightingPotentialGrids
void SetWeightingPotentialGrids(const double xmin, const double xmax, const double xsteps, const double ymin, const double ymax, const double ysteps, const double zmin, const double zmax, const double zsteps, const double tmin, const double tmax, const double tsteps)
Definition: ComponentParallelPlate.cc:718
Garfield::ComponentParallelPlate::GetVoltageRange
bool GetVoltageRange(double &vmin, double &vmax) override
Calculate the voltage range [V].
Definition: ComponentParallelPlate.cc:464
Garfield::ComponentParallelPlate::GetBoundingBox
bool GetBoundingBox(double &xmin, double &ymin, double &zmin, double &xmax, double &ymax, double &zmax) override
Get the bounding box coordinates.
Definition: ComponentParallelPlate.cc:40
Garfield::ComponentParallelPlate::AddPixel
void AddPixel(double x, double y, double lx, double ly, const std::string &label)
Definition: ComponentParallelPlate.cc:586
Garfield::Component
Abstract base class for components.
Definition: Component.hh:13
Garfield::Component::m_debug
bool m_debug
Switch on/off debugging messages.
Definition: Component.hh:341
Garfield::Component::m_className
std::string m_className
Class name.
Definition: Component.hh:329
Garfield::Component::m_geometry
Geometry * m_geometry
Pointer to the geometry.
Definition: Component.hh:332
Garfield::Geometry::GetBoundingBox
virtual bool GetBoundingBox(double &xmin, double &ymin, double &zmin, double &xmax, double &ymax, double &zmax)=0
Get the bounding box (envelope of the geometry).
Garfield::Geometry::GetMedium
virtual Medium * GetMedium(const double x, const double y, const double z, const bool tesselated=false) const =0
Retrieve the medium at a given point.
Garfield::Medium
Abstract base class for media.
Definition: Medium.hh:13
Heed::cos
DoubleAc cos(const DoubleAc &f)
Definition: DoubleAc.cpp:432
Heed::exp
DoubleAc exp(const DoubleAc &f)
Definition: DoubleAc.cpp:377
Heed::sin
DoubleAc sin(const DoubleAc &f)
Definition: DoubleAc.cpp:384