Garfield++ 5.0
A toolkit for the detailed simulation of particle detectors based on ionisation measurement in gases and semiconductors
Loading...
Searching...
No Matches
ViewBase.cc
Go to the documentation of this file.
1#include <iostream>
2#include <cstdio>
3#include <cmath>
4#include <limits>
5
6#include <TROOT.h>
7#include <TGraph.h>
8
9#include "Garfield/Sensor.hh"
10#include "Garfield/Component.hh"
11#include "Garfield/Plotting.hh"
12#include "Garfield/GarfieldConstants.hh"
13#include "Garfield/ViewBase.hh"
14
15namespace {
16
17bool Invert(std::array<std::array<double, 3>, 3>& a) {
18
19 // Compute cofactors.
20 const double c11 = a[1][1] * a[2][2] - a[1][2] * a[2][1];
21 const double c12 = a[1][2] * a[2][0] - a[1][0] * a[2][2];
22 const double c13 = a[1][0] * a[2][1] - a[1][1] * a[2][0];
23 const double c21 = a[2][1] * a[0][2] - a[2][2] * a[0][1];
24 const double c22 = a[2][2] * a[0][0] - a[2][0] * a[0][2];
25 const double c23 = a[2][0] * a[0][1] - a[2][1] * a[0][0];
26 const double c31 = a[0][1] * a[1][2] - a[0][2] * a[1][1];
27 const double c32 = a[0][2] * a[1][0] - a[0][0] * a[1][2];
28 const double c33 = a[0][0] * a[1][1] - a[0][1] * a[1][0];
29 const double t1 = fabs(a[0][0]);
30 const double t2 = fabs(a[1][0]);
31 const double t3 = fabs(a[2][0]);
32 double det = 0.;
33 double pivot = 0.;
34 if (t2 < t1 && t3 < t1) {
35 pivot = a[0][0];
36 det = c22 * c33 - c23 * c32;
37 } else if (t1 < t2 && t3 < t2) {
38 pivot = a[1][0];
39 det = c13 * c32 - c12 * c33;
40 } else {
41 pivot = a[2][0];
42 det = c23 * c12 - c22 * c13;
43 }
44 if (det == 0.) return false;
45 const double s = pivot / det;
46 a[0][0] = s * c11;
47 a[0][1] = s * c21;
48 a[0][2] = s * c31;
49 a[1][0] = s * c12;
50 a[1][1] = s * c22;
51 a[1][2] = s * c32;
52 a[2][0] = s * c13;
53 a[2][1] = s * c23;
54 a[2][2] = s * c33;
55 return true;
56}
57
58std::string Fmt(const double x) {
59 char buf[100];
60 snprintf(buf, 100, "%g", x);
61 return std::string(buf);
62}
63
64}
65
66namespace Garfield {
67
68ViewBase::ViewBase(const std::string& name) :
69 m_className(name) {
70
71 plottingEngine.SetDefaultStyle();
72}
73
74TPad* ViewBase::GetCanvas() {
75 if (!m_pad) {
76 std::string name = FindUnusedCanvasName("c" + m_className);
77 if (!m_canvas) m_canvas.reset(new TCanvas(name.c_str(), ""));
78 m_pad = m_canvas.get();
79 }
80 return m_pad;
81}
82
83bool ViewBase::RangeSet(TVirtualPad* pad) {
84 if (!pad) return false;
85 if (pad->GetListOfPrimitives()->GetSize() == 0 &&
86 pad->GetX1() == 0 && pad->GetX2() == 1 &&
87 pad->GetY1() == 0 && pad->GetY2() == 1) {
88 return false;
89 }
90 return true;
91}
92
93void ViewBase::SetRange(TVirtualPad* pad, const double x0, const double y0,
94 const double x1, const double y1) {
95 if (!pad) return;
96 const double bm = pad->GetBottomMargin();
97 const double lm = pad->GetLeftMargin();
98 const double rm = pad->GetRightMargin();
99 const double tm = pad->GetTopMargin();
100 const double dx = x1 - x0;
101 const double dy = y1 - y0;
102 pad->Range(x0 - dx * (lm / (1. - rm - lm)),
103 y0 - dy * (bm / (1. - tm - lm)),
104 x1 + dx * (rm / (1. - rm - lm)),
105 y1 + dy * (tm / (1. - tm - lm)));
106
107}
108
109void ViewBase::SetArea(const double xmin, const double ymin,
110 const double xmax, const double ymax) {
111 // Check range, assign if non-null.
112 if (xmin == xmax || ymin == ymax) {
113 std::cerr << m_className << "::SetArea: Null area is not permitted.\n"
114 << " " << xmin << " < x < " << xmax << "\n"
115 << " " << ymin << " < y < " << ymax << "\n";
116 return;
117 }
118 m_xMinPlot = std::min(xmin, xmax);
119 m_yMinPlot = std::min(ymin, ymax);
120 m_xMaxPlot = std::max(xmin, xmax);
121 m_yMaxPlot = std::max(ymin, ymax);
122 m_userPlotLimits = true;
123}
124
125void ViewBase::SetArea(const double xmin, const double ymin, const double zmin,
126 const double xmax, const double ymax,
127 const double zmax) {
128 // Check range, assign if non-null
129 if (xmin == xmax || ymin == ymax || zmin == zmax) {
130 std::cerr << m_className << "::SetArea: Null area range not permitted.\n";
131 return;
132 }
133 m_xMinBox = std::min(xmin, xmax);
134 m_yMinBox = std::min(ymin, ymax);
135 m_zMinBox = std::min(zmin, zmax);
136 m_xMaxBox = std::max(xmin, xmax);
137 m_yMaxBox = std::max(ymin, ymax);
138 m_zMaxBox = std::max(zmin, zmax);
139 m_userBox = true;
140}
141
142void ViewBase::SetPlane(const double fx, const double fy, const double fz,
143 const double x0, const double y0, const double z0) {
144 // Calculate two in-plane vectors for the normal vector
145 const double fnorm = sqrt(fx * fx + fy * fy + fz * fz);
146 if (fnorm > 0 && fx * fx + fz * fz > 0) {
147 const double fxz = sqrt(fx * fx + fz * fz);
148 m_proj[0][0] = fz / fxz;
149 m_proj[0][1] = 0;
150 m_proj[0][2] = -fx / fxz;
151 m_proj[1][0] = -fx * fy / (fxz * fnorm);
152 m_proj[1][1] = (fx * fx + fz * fz) / (fxz * fnorm);
153 m_proj[1][2] = -fy * fz / (fxz * fnorm);
154 m_proj[2][0] = x0;
155 m_proj[2][1] = y0;
156 m_proj[2][2] = z0;
157 } else if (fnorm > 0 && fy * fy + fz * fz > 0) {
158 const double fyz = sqrt(fy * fy + fz * fz);
159 m_proj[0][0] = (fy * fy + fz * fz) / (fyz * fnorm);
160 m_proj[0][1] = -fx * fz / (fyz * fnorm);
161 m_proj[0][2] = -fy * fz / (fyz * fnorm);
162 m_proj[1][0] = 0;
163 m_proj[1][1] = fz / fyz;
164 m_proj[1][2] = -fy / fyz;
165 m_proj[2][0] = x0;
166 m_proj[2][1] = y0;
167 m_proj[2][2] = z0;
168 } else {
169 std::cout << m_className << "::SetPlane:\n"
170 << " Normal vector has zero norm. No new projection set.\n";
171 }
172
173 // Store the plane description
174 m_plane[0] = fx;
175 m_plane[1] = fy;
176 m_plane[2] = fz;
177 m_plane[3] = fx * x0 + fy * y0 + fz * z0;
178
179 UpdateProjectionMatrix();
180}
181
182void ViewBase::SetPlane(const double fx, const double fy, const double fz,
183 const double x0, const double y0, const double z0,
184 const double hx, const double hy, const double hz) {
185
186 const double fnorm = sqrt(fx * fx + fy * fy + fz * fz);
187 if (fnorm < Small) {
188 std::cout << m_className << "::SetPlane:\n"
189 << " Normal vector has zero norm. No new projection set.\n";
190 return;
191 }
192 // Normalise the vector.
193 const double wx = fx / fnorm;
194 const double wy = fy / fnorm;
195 const double wz = fz / fnorm;
196 // Store the plane description.
197 m_plane[0] = wx;
198 m_plane[1] = wy;
199 m_plane[2] = wz;
200 m_plane[3] = wx * x0 + wy * y0 + wz * z0;
201
202 double d = hx * wx + hy * wy + hz * wz;
203 double ux = hx - d * wx;
204 double uy = hy - d * wy;
205 double uz = hz - d * wz;
206 double unorm = std::sqrt(ux * ux + uy * uy + uz * uz);
207 if (unorm < 1.e-10) {
208 // Wrong in-plane x hint (close to norm).
209 if (fy * fy + fz * fz > 0) {
210 // Taking global x as in-plane x hint.
211 ux = 1;
212 uy = 0;
213 uz = 0;
214 } else {
215 // Taking global y as in-plane x hint.
216 ux = 0;
217 uy = 1;
218 uz = 0;
219 }
220 d = ux * wx + uy * wy + uz * wz;
221 ux -= d * wx;
222 uy -= d * wy;
223 uz -= d * wz;
224 unorm = std::sqrt(ux * ux + uy * uy + uz * uz);
225 }
226 ux /= unorm;
227 uy /= unorm;
228 uz /= unorm;
229
230 m_prmat[0][0] = ux;
231 m_prmat[1][0] = uy;
232 m_prmat[2][0] = uz;
233 // In-plane y = cross product [z,x]
234 m_prmat[0][1] = wy * uz - wz * uy;
235 m_prmat[1][1] = wz * ux - wx * uz;
236 m_prmat[2][1] = wx * uy - wy * ux;
237 m_prmat[0][2] = wx;
238 m_prmat[1][2] = wy;
239 m_prmat[2][2] = wz;
240
241 for (unsigned int i = 0; i < 3; ++i) {
242 m_proj[0][i] = m_prmat[i][0];
243 m_proj[1][i] = m_prmat[i][1];
244 }
245 m_proj[2][0] = x0;
246 m_proj[2][1] = y0;
247 m_proj[2][2] = z0;
248 if (!Invert(m_prmat)) {
249 std::cerr << m_className << "::SetPlane:\n"
250 << " Inversion failed; reset to default.\n";
251 SetPlaneXY();
252 }
253 if (m_debug) {
254 std::cout << m_className << "::SetPlane:\n PRMAT:\n";
255 for (size_t i = 0; i < 3; ++i) {
256 std::printf(" %10.5f %10.5f %10.5f\n",
257 m_prmat[i][0], m_prmat[i][1], m_prmat[i][2]);
258 }
259 std::cout << " PROJ:\n";
260 for (size_t i = 0; i < 3; ++i) {
261 std::printf(" %10.5f %10.5f %10.5f\n",
262 m_proj[i][0], m_proj[i][1], m_proj[i][2]);
263 }
264 std::cout << " PLANE:\n";
265 std::printf(" %10.5f %10.5f %10.5f %10.5f\n",
266 m_plane[0], m_plane[1], m_plane[2], m_plane[3]);
267 }
268}
269
270void ViewBase::Rotate(const double theta) {
271 // Rotate the axes
272 double auxu[3], auxv[3];
273 const double ctheta = cos(theta);
274 const double stheta = sin(theta);
275 for (int i = 0; i < 3; ++i) {
276 auxu[i] = ctheta * m_proj[0][i] - stheta * m_proj[1][i];
277 auxv[i] = stheta * m_proj[0][i] + ctheta * m_proj[1][i];
278 }
279 for (int i = 0; i < 3; ++i) {
280 m_proj[0][i] = auxu[i];
281 m_proj[1][i] = auxv[i];
282 }
283
284 UpdateProjectionMatrix();
285}
286
287void ViewBase::SetPlaneXY() {
288 m_proj = {{{1, 0, 0}, {0, 1, 0}, {0, 0, 0}}};
289 m_plane = {0, 0, 1, 0};
290 m_prmat = {{{1, 0, 0}, {0, 1, 0}, {0, 0, 1}}};
291}
292
293void ViewBase::SetPlaneXZ() {
294 m_proj = {{{1, 0, 0}, {0, 0, 1}, {0, 0, 0}}};
295 m_plane = {0, 1, 0, 0};
296 UpdateProjectionMatrix();
297}
298
299void ViewBase::SetPlaneYZ() {
300 m_proj = {{{0, 1, 0}, {0, 0, 1}, {0, 0, 0}}};
301 m_plane = {1, 0, 0, 0};
302 UpdateProjectionMatrix();
303}
304
305void ViewBase::SetPlaneZX() {
306 m_proj = {{{0, 0, 1}, {1, 0, 0}, {0, 0, 0}}};
307 m_plane = {0, 1, 0, 0};
308 UpdateProjectionMatrix();
309}
310
311void ViewBase::SetPlaneZY() {
312 m_proj = {{{0, 0, 1}, {0, 1, 0}, {0, 0, 0}}};
313 m_plane = {1, 0, 0, 0};
314 UpdateProjectionMatrix();
315}
316
317std::string ViewBase::FindUnusedFunctionName(const std::string& s) {
318 int idx = 0;
319 std::string fname = s + "_0";
320 while (gROOT->GetListOfFunctions()->FindObject(fname.c_str())) {
321 ++idx;
322 fname = s + "_" + std::to_string(idx);
323 }
324 return fname;
325}
326
327std::string ViewBase::FindUnusedHistogramName(const std::string& s) {
328 int idx = 0;
329 std::string hname = s + "_0";
330 while (gDirectory->GetList()->FindObject(hname.c_str())) {
331 ++idx;
332 hname = s + "_" + std::to_string(idx);
333 }
334 return hname;
335}
336
337std::string ViewBase::FindUnusedCanvasName(const std::string& s) {
338 int idx = 0;
339 std::string hname = s + "_0";
340 while (gROOT->GetListOfCanvases()->FindObject(hname.c_str())) {
341 ++idx;
342 hname = s + "_" + std::to_string(idx);
343 }
344 return hname;
345}
346
347void ViewBase::UpdateProjectionMatrix() {
348
349 m_prmat[0][0] = m_proj[0][0];
350 m_prmat[1][0] = m_proj[0][1];
351 m_prmat[2][0] = m_proj[0][2];
352 m_prmat[0][1] = m_proj[1][0];
353 m_prmat[1][1] = m_proj[1][1];
354 m_prmat[2][1] = m_proj[1][2];
355 const double vnorm = sqrt(m_plane[0] * m_plane[0] +
356 m_plane[1] * m_plane[1] +
357 m_plane[2] * m_plane[2]);
358 if (vnorm <= 0.) {
359 std::cerr << m_className << "::UpdateProjectionMatrix:\n"
360 << " Zero norm vector; reset to default.\n";
361 SetPlaneXY();
362 return;
363 }
364 m_prmat[0][2] = m_plane[0] / vnorm;
365 m_prmat[1][2] = m_plane[1] / vnorm;
366 m_prmat[2][2] = m_plane[2] / vnorm;
367 if (!Invert(m_prmat)) {
368 std::cerr << m_className << "::UpdateProjectionMatrix:\n"
369 << " Inversion failed; reset to default.\n";
370 SetPlaneXY();
371 }
372}
373
374void ViewBase::Clip(const std::array<float, 3>& x0,
375 const std::array<float, 3>& x1,
376 std::array<float, 3>& xc) const {
377
378 xc.fill(0.);
379 const bool in0 = InBox(x0);
380 const bool in1 = InBox(x1);
381 if (in0 == in1) return;
382 xc = in0 ? x1 : x0;
383 const std::array<float, 3> dx = {x1[0] - x0[0], x1[1] - x0[1],
384 x1[2] - x0[2]};
385 std::array<double, 3> bmin = {m_xMinBox, m_yMinBox, m_zMinBox};
386 std::array<double, 3> bmax = {m_xMaxBox, m_yMaxBox, m_zMaxBox};
387 for (size_t i = 0; i < 3; ++i) {
388 if (dx[i] == 0. || (xc[i] >= bmin[i] && xc[i] <= bmax[i])) continue;
389 const double b = xc[i] < bmin[i] ? bmin[i] : bmax[i];
390 const double s = (b - xc[i]) / dx[i];
391 xc[i] = b;
392 for (size_t j = 0; j < 3; ++j) {
393 if (j != i) xc[j] += dx[j] * s;
394 }
395 }
396}
397
398void ViewBase::DrawLine(const std::vector<std::array<float, 3> >& xl,
399 const short col, const short lw) {
400
401 const size_t nP = xl.size();
402 if (nP < 2) return;
403
404 TGraph gr;
405 gr.SetLineColor(col);
406 gr.SetLineWidth(lw);
407
408 std::vector<float> xgr;
409 std::vector<float> ygr;
410 auto x0 = xl[0];
411 bool in0 = InBox(x0);
412 if (in0) {
413 float xp = 0., yp = 0.;
414 ToPlane(x0[0], x0[1], x0[2], xp, yp);
415 xgr.push_back(xp);
416 ygr.push_back(yp);
417 }
418 for (unsigned int j = 1; j < nP; ++j) {
419 auto x1 = xl[j];
420 bool in1 = InBox(x1);
421 if (in1 != in0) {
422 float xp = 0., yp = 0.;
423 std::array<float, 3> xc;
424 Clip(x0, x1, xc);
425 ToPlane(xc[0], xc[1], xc[2], xp, yp);
426 xgr.push_back(xp);
427 ygr.push_back(yp);
428 }
429 if (in1) {
430 float xp = 0., yp = 0.;
431 ToPlane(x1[0], x1[1], x1[2], xp, yp);
432 xgr.push_back(xp);
433 ygr.push_back(yp);
434 } else if (!xgr.empty()) {
435 gr.DrawGraph(xgr.size(), xgr.data(), ygr.data(), "Lsame");
436 xgr.clear();
437 ygr.clear();
438 }
439 x0 = x1;
440 in0 = in1;
441 }
442 if (!xgr.empty()) {
443 gr.DrawGraph(xgr.size(), xgr.data(), ygr.data(), "Lsame");
444 }
445}
446
447std::string ViewBase::LabelX() {
448
449 std::string xLabel = "";
450 constexpr double tol = 1.e-4;
451 // x portion
452 if (fabs(m_proj[0][0] - 1) < tol) {
453 xLabel = "#it{x}";
454 } else if (fabs(m_proj[0][0] + 1) < tol) {
455 xLabel = "#minus#it{x}";
456 } else if (fabs(m_proj[0][0]) > tol) {
457 xLabel = Fmt(m_proj[0][0]) + " #it{x}";
458 }
459
460 // y portion
461 if (!xLabel.empty()) {
462 if (m_proj[0][1] < -tol) {
463 xLabel += " #minus ";
464 } else if (m_proj[0][1] > tol) {
465 xLabel += " #plus ";
466 }
467 if (fabs(m_proj[0][1] - 1) < tol || fabs(m_proj[0][1] + 1) < tol) {
468 xLabel += "#it{y}";
469 } else if (fabs(m_proj[0][1]) > tol) {
470 xLabel += Fmt(fabs(m_proj[0][1])) + " #it{y}";
471 }
472 } else {
473 if (fabs(m_proj[0][1] - 1) < tol) {
474 xLabel = "#it{y}";
475 } else if (fabs(m_proj[0][1] + 1) < tol) {
476 xLabel = "#minus#it{y}";
477 } else if (fabs(m_proj[0][1]) > tol) {
478 xLabel = Fmt(m_proj[0][1]) + " #it{y}";
479 }
480 }
481
482 // z portion
483 if (!xLabel.empty()) {
484 if (m_proj[0][2] < -tol) {
485 xLabel += " #minus ";
486 } else if (m_proj[0][2] > tol) {
487 xLabel += " #plus ";
488 }
489 if (fabs(m_proj[0][2] - 1) < tol || fabs(m_proj[0][2] + 1) < tol) {
490 xLabel += "#it{z}";
491 } else if (fabs(m_proj[0][2]) > tol) {
492 xLabel += Fmt(fabs(m_proj[0][2])) + " #it{z}";
493 }
494 } else {
495 if (fabs(m_proj[0][2] - 1) < tol) {
496 xLabel = "#it{z}";
497 } else if (fabs(m_proj[0][2] + 1) < tol) {
498 xLabel = "#minus#it{z}";
499 } else if (fabs(m_proj[0][2]) > tol) {
500 xLabel = Fmt(m_proj[0][2]) + " #it{z}";
501 }
502 }
503
504 // Unit
505 xLabel += " [cm]";
506 return xLabel;
507
508}
509
510std::string ViewBase::LabelY() {
511
512 std::string yLabel = "";
513 constexpr double tol = 1.e-4;
514 // x portion
515 if (fabs(m_proj[1][0] - 1) < tol) {
516 yLabel = "#it{x}";
517 } else if (fabs(m_proj[1][0] + 1) < tol) {
518 yLabel = "#minus#it{x}";
519 } else if (fabs(m_proj[1][0]) > tol) {
520 yLabel = Fmt(m_proj[1][0]) + " #it{x}";
521 }
522
523 // y portion
524 if (!yLabel.empty()) {
525 if (m_proj[1][1] < -tol) {
526 yLabel += " #minus ";
527 } else if (m_proj[1][1] > tol) {
528 yLabel += " #plus ";
529 }
530 if (fabs(m_proj[1][1] - 1) < tol || fabs(m_proj[1][1] + 1) < tol) {
531 yLabel += "#it{y}";
532 } else if (fabs(m_proj[1][1]) > tol) {
533 yLabel += Fmt(fabs(m_proj[1][1])) + " #it{y}";
534 }
535 } else {
536 if (fabs(m_proj[1][1] - 1) < tol) {
537 yLabel = "#it{y}";
538 } else if (fabs(m_proj[1][1] + 1) < tol) {
539 yLabel = "#minus#it{y}";
540 } else if (fabs(m_proj[1][1]) > tol) {
541 yLabel = Fmt(m_proj[1][1]) + " #it{y}";
542 }
543 }
544
545 // z portion
546 if (!yLabel.empty()) {
547 if (m_proj[1][2] < -tol) {
548 yLabel += " #minus ";
549 } else if (m_proj[1][2] > tol) {
550 yLabel += " #plus ";
551 }
552 if (fabs(m_proj[1][2] - 1) < tol || fabs(m_proj[1][2] + 1) < tol) {
553 yLabel += "#it{z}";
554 } else if (fabs(m_proj[1][2]) > tol) {
555 yLabel += Fmt(fabs(m_proj[1][2])) + " #it{z}";
556 }
557 } else {
558 if (fabs(m_proj[1][2] - 1) < tol) {
559 yLabel = "#it{z}";
560 } else if (fabs(m_proj[1][2] + 1) < tol) {
561 yLabel = "#minus#it{z}";
562 } else if (fabs(m_proj[1][2]) > tol) {
563 yLabel = Fmt(m_proj[1][2]) + " #it{z}";
564 }
565 }
566
567 // Unit
568 yLabel += " [cm]";
569 return yLabel;
570}
571
572std::string ViewBase::PlaneDescription() {
573
574 std::string description;
575
576 constexpr double tol = 1.e-4;
577 // x portion
578 if (fabs(m_plane[0] - 1) < tol) {
579 description = "x";
580 } else if (fabs(m_plane[0] + 1) < tol) {
581 description = "-x";
582 } else if (fabs(m_plane[0]) > tol) {
583 description = Fmt(m_plane[0]) + " x";
584 }
585
586 // y portion
587 if (!description.empty()) {
588 if (m_plane[1] < -tol) {
589 description += " - ";
590 } else if (m_plane[1] > tol) {
591 description += " + ";
592 }
593 if (fabs(m_plane[1] - 1) < tol || fabs(m_plane[1] + 1) < tol) {
594 description += "y";
595 } else if (fabs(m_plane[1]) > tol) {
596 description += Fmt(fabs(m_plane[1])) + " y";
597 }
598 } else {
599 if (fabs(m_plane[1] - 1) < tol) {
600 description = "y";
601 } else if (fabs(m_plane[1] + 1) < tol) {
602 description = "-y";
603 } else if (fabs(m_plane[1]) > tol) {
604 description = Fmt(m_plane[1]) + " y";
605 }
606 }
607
608 // z portion
609 if (!description.empty()) {
610 if (m_plane[2] < -tol) {
611 description += " - ";
612 } else if (m_plane[2] > tol) {
613 description += " + ";
614 }
615 if (fabs(m_plane[2] - 1) < tol || fabs(m_plane[2] + 1) < tol) {
616 description += "z";
617 } else if (fabs(m_plane[2]) > tol) {
618 description += Fmt(fabs(m_plane[2])) + " z";
619 }
620 } else {
621 if (fabs(m_plane[2] - 1) < tol) {
622 description = "z";
623 } else if (fabs(m_plane[2] + 1) < tol) {
624 description = "-z";
625 } else if (fabs(m_plane[2]) > tol) {
626 description = Fmt(m_plane[2]) + " z";
627 }
628 }
629
630 // Constant
631 description += " = " + Fmt(m_plane[3]);
632 return description;
633}
634
635bool ViewBase::PlotLimits(Sensor* sensor,
636 double& xmin, double& ymin,
637 double& xmax, double& ymax) const {
638
639 if (!sensor) return false;
640 // Try to get the area/bounding box from the sensor/component.
641 std::array<double, 3> bbmin;
642 std::array<double, 3> bbmax;
643 if (!sensor->GetArea(bbmin[0], bbmin[1], bbmin[2],
644 bbmax[0], bbmax[1], bbmax[2])) {
645 std::cerr << m_className << "::PlotLimits:\n"
646 << " Sensor area is not defined.\n"
647 << " Please set the plot limits explicitly (SetArea).\n";
648 return false;
649 }
650 return PlotLimits(bbmin, bbmax, xmin, ymin, xmax, ymax);
651}
652
653bool ViewBase::PlotLimits(Component* cmp,
654 double& xmin, double& ymin,
655 double& xmax, double& ymax) const {
656
657 if (!cmp) return false;
658 // Try to get the area/bounding box from the sensor/component.
659 std::array<double, 3> bbmin;
660 std::array<double, 3> bbmax;
661 if (!cmp->GetBoundingBox(bbmin[0], bbmin[1], bbmin[2],
662 bbmax[0], bbmax[1], bbmax[2])) {
663 std::cerr << m_className << "::PlotLimits:\n"
664 << " Bounding box of the component is not defined.\n"
665 << " Please set the plot limits explicitly (SetArea).\n";
666 return false;
667 }
668 if (std::isinf(bbmin[0]) || std::isinf(bbmax[0]) ||
669 std::isinf(bbmin[1]) || std::isinf(bbmax[1]) ||
670 std::isinf(bbmin[2]) || std::isinf(bbmax[2])) {
671 std::array<double, 3> cellmin = {0., 0., 0.};
672 std::array<double, 3> cellmax = {0., 0., 0.};
673 if (!cmp->GetElementaryCell(cellmin[0], cellmin[1], cellmin[2],
674 cellmax[0], cellmax[1], cellmax[2])) {
675 std::cerr << m_className << "::PlotLimits:\n"
676 << " Cell boundaries are not defined.\n"
677 << " Please set the plot limits explicitly (SetArea).\n";
678 return false;
679 }
680 for (size_t i = 0; i < 3; ++i) {
681 if (std::isinf(bbmin[i]) || std::isinf(bbmax[i])) {
682 bbmin[i] = cellmin[i];
683 bbmax[i] = cellmax[i];
684 }
685 }
686 }
687 return PlotLimits(bbmin, bbmax, xmin, ymin, xmax, ymax);
688}
689
690bool ViewBase::PlotLimitsFromUserBox(double& xmin, double& ymin,
691 double& xmax, double& ymax) const {
692
693 std::array<double, 3> bbmin = {m_xMinBox, m_yMinBox, m_zMinBox};
694 std::array<double, 3> bbmax = {m_xMaxBox, m_yMaxBox, m_zMaxBox};
695 return PlotLimits(bbmin, bbmax, xmin, ymin, xmax, ymax);
696}
697
698bool ViewBase::PlotLimits(std::array<double, 3>& bbmin,
699 std::array<double, 3>& bbmax,
700 double& xmin, double& ymin,
701 double& xmax, double& ymax) const {
702 constexpr double tol = 1.e-4;
703 double umin[2] = {-std::numeric_limits<double>::max(),
704 -std::numeric_limits<double>::max()};
705 double umax[2] = {std::numeric_limits<double>::max(),
706 std::numeric_limits<double>::max()};
707 for (unsigned int i = 0; i < 3; ++i) {
708 bbmin[i] -= m_proj[2][i];
709 bbmax[i] -= m_proj[2][i];
710 for (unsigned int j = 0; j < 2; ++j) {
711 if (fabs(m_proj[j][i]) < tol) continue;
712 const double t1 = bbmin[i] / m_proj[j][i];
713 const double t2 = bbmax[i] / m_proj[j][i];
714 const double tmin = std::min(t1, t2);
715 const double tmax = std::max(t1, t2);
716 if (tmin > umin[j] && tmin < umax[j]) umin[j] = tmin;
717 if (tmax < umax[j] && tmax > umin[j]) umax[j] = tmax;
718 }
719 }
720 xmin = umin[0];
721 xmax = umax[0];
722 ymin = umin[1];
723 ymax = umax[1];
724 return true;
725}
726
727}
Garfield::Component
Abstract base class for components.
Definition Component.hh:13
Garfield::Component::GetBoundingBox
virtual bool GetBoundingBox(double &xmin, double &ymin, double &zmin, double &xmax, double &ymax, double &zmax)
Get the bounding box coordinates.
Definition Component.cc:120
Garfield::Component::GetElementaryCell
virtual bool GetElementaryCell(double &xmin, double &ymin, double &zmin, double &xmax, double &ymax, double &zmax)
Get the coordinates of the elementary cell.
Definition Component.cc:126
Garfield::Sensor::GetArea
bool GetArea(double &xmin, double &ymin, double &zmin, double &xmax, double &ymax, double &zmax)
Return the current user area.
Definition Sensor.cc:234
Garfield::ViewBase::m_plane
std::array< double, 4 > m_plane
Definition ViewBase.hh:103
Garfield::ViewBase::LabelY
std::string LabelY()
Definition ViewBase.cc:510
Garfield::ViewBase::SetRange
static void SetRange(TVirtualPad *pad, const double x0, const double y0, const double x1, const double y1)
Definition ViewBase.cc:93
Garfield::ViewBase::LabelX
std::string LabelX()
Definition ViewBase.cc:447
Garfield::ViewBase::m_prmat
std::array< std::array< double, 3 >, 3 > m_prmat
Definition ViewBase.hh:105
Garfield::ViewBase::Clip
void Clip(const std::array< float, 3 > &x0, const std::array< float, 3 > &x1, std::array< float, 3 > &xc) const
Definition ViewBase.cc:374
Garfield::ViewBase::InBox
bool InBox(const std::array< T, 3 > &x) const
Definition ViewBase.hh:119
Garfield::ViewBase::DrawLine
void DrawLine(const std::vector< std::array< float, 3 > > &xl, const short col, const short lw)
Definition ViewBase.cc:398
Garfield::ViewBase::FindUnusedHistogramName
static std::string FindUnusedHistogramName(const std::string &s)
Find an unused histogram name.
Definition ViewBase.cc:327
Garfield::ViewBase::m_className
std::string m_className
Definition ViewBase.hh:82
Garfield::ViewBase::SetPlaneZX
void SetPlaneZX()
Set the viewing plane to z-x.
Definition ViewBase.cc:305
Garfield::ViewBase::SetPlaneYZ
void SetPlaneYZ()
Set the viewing plane to y-z.
Definition ViewBase.cc:299
Garfield::ViewBase::FindUnusedCanvasName
static std::string FindUnusedCanvasName(const std::string &s)
Find an unused canvas name.
Definition ViewBase.cc:337
Garfield::ViewBase::m_proj
std::array< std::array< double, 3 >, 3 > m_proj
Definition ViewBase.hh:100
Garfield::ViewBase::GetCanvas
TPad * GetCanvas()
Retrieve the canvas.
Definition ViewBase.cc:74
Garfield::ViewBase::RangeSet
static bool RangeSet(TVirtualPad *)
Definition ViewBase.cc:83
Garfield::ViewBase::SetPlaneXZ
void SetPlaneXZ()
Set the viewing plane to x-z.
Definition ViewBase.cc:293
Garfield::ViewBase::PlotLimits
bool PlotLimits(Sensor *sensor, double &xmin, double &ymin, double &xmax, double &ymax) const
Definition ViewBase.cc:635
Garfield::ViewBase::SetPlaneXY
void SetPlaneXY()
Set the viewing plane to x-y.
Definition ViewBase.cc:287
Garfield::ViewBase::ToPlane
void ToPlane(const T x, const T y, const T z, T &xp, T &yp) const
Definition ViewBase.hh:113
Garfield::ViewBase::PlaneDescription
std::string PlaneDescription()
Definition ViewBase.cc:572
Garfield::ViewBase::SetPlaneZY
void SetPlaneZY()
Set the viewing plane to z-y.
Definition ViewBase.cc:311
Garfield::ViewBase::PlotLimitsFromUserBox
bool PlotLimitsFromUserBox(double &xmin, double &ymin, double &xmax, double &ymax) const
Definition ViewBase.cc:690
Garfield::ViewBase::FindUnusedFunctionName
static std::string FindUnusedFunctionName(const std::string &s)
Find an unused function name.
Definition ViewBase.cc:317
Garfield::ViewBase::Rotate
void Rotate(const double angle)
Rotate the viewing plane (angle in radian).
Definition ViewBase.cc:270
Garfield::ViewBase::ViewBase
ViewBase()=delete
Default constructor.
Garfield::ViewBase::UpdateProjectionMatrix
void UpdateProjectionMatrix()
Definition ViewBase.cc:347
Garfield::ViewBase::SetPlane
virtual void SetPlane(const double fx, const double fy, const double fz, const double x0, const double y0, const double z0)
Definition ViewBase.cc:142
Garfield::plottingEngine
PlottingEngine plottingEngine
Definition PlottingEngine.cc:13
Heed::fabs
DoubleAc fabs(const DoubleAc &f)
Definition DoubleAc.h:615