Garfield++ 4.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 sprintf(buf, "%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(TPad* pad) {
84
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(TPad* 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}
254
255void ViewBase::Rotate(const double theta) {
256 // Rotate the axes
257 double auxu[3], auxv[3];
258 const double ctheta = cos(theta);
259 const double stheta = sin(theta);
260 for (int i = 0; i < 3; ++i) {
261 auxu[i] = ctheta * m_proj[0][i] - stheta * m_proj[1][i];
262 auxv[i] = stheta * m_proj[0][i] + ctheta * m_proj[1][i];
263 }
264 for (int i = 0; i < 3; ++i) {
265 m_proj[0][i] = auxu[i];
266 m_proj[1][i] = auxv[i];
267 }
268
269 UpdateProjectionMatrix();
270}
271
272void ViewBase::SetPlaneXY() {
273 m_proj = {{{1, 0, 0}, {0, 1, 0}, {0, 0, 0}}};
274 m_plane = {0, 0, 1, 0};
275 m_prmat = {{{1, 0, 0}, {0, 1, 0}, {0, 0, 1}}};
276}
277
278void ViewBase::SetPlaneXZ() {
279 m_proj = {{{1, 0, 0}, {0, 0, 1}, {0, 0, 0}}};
280 m_plane = {0, 1, 0, 0};
281 UpdateProjectionMatrix();
282}
283
284void ViewBase::SetPlaneYZ() {
285 m_proj = {{{0, 1, 0}, {0, 0, 1}, {0, 0, 0}}};
286 m_plane = {1, 0, 0, 0};
287 UpdateProjectionMatrix();
288}
289
290std::string ViewBase::FindUnusedFunctionName(const std::string& s) {
291 int idx = 0;
292 std::string fname = s + "_0";
293 while (gROOT->GetListOfFunctions()->FindObject(fname.c_str())) {
294 ++idx;
295 fname = s + "_" + std::to_string(idx);
296 }
297 return fname;
298}
299
300std::string ViewBase::FindUnusedHistogramName(const std::string& s) {
301 int idx = 0;
302 std::string hname = s + "_0";
303 while (gDirectory->GetList()->FindObject(hname.c_str())) {
304 ++idx;
305 hname = s + "_" + std::to_string(idx);
306 }
307 return hname;
308}
309
310std::string ViewBase::FindUnusedCanvasName(const std::string& s) {
311 int idx = 0;
312 std::string hname = s + "_0";
313 while (gROOT->GetListOfCanvases()->FindObject(hname.c_str())) {
314 ++idx;
315 hname = s + "_" + std::to_string(idx);
316 }
317 return hname;
318}
319
320void ViewBase::UpdateProjectionMatrix() {
321
322 m_prmat[0][0] = m_proj[0][0];
323 m_prmat[1][0] = m_proj[0][1];
324 m_prmat[2][0] = m_proj[0][2];
325 m_prmat[0][1] = m_proj[1][0];
326 m_prmat[1][1] = m_proj[1][1];
327 m_prmat[2][1] = m_proj[1][2];
328 const double vnorm = sqrt(m_plane[0] * m_plane[0] +
329 m_plane[1] * m_plane[1] +
330 m_plane[2] * m_plane[2]);
331 if (vnorm <= 0.) {
332 std::cerr << m_className << "::UpdateProjectionMatrix:\n"
333 << " Zero norm vector; reset to default.\n";
334 SetPlaneXY();
335 return;
336 }
337 m_prmat[0][2] = m_plane[0] / vnorm;
338 m_prmat[1][2] = m_plane[1] / vnorm;
339 m_prmat[2][2] = m_plane[2] / vnorm;
340 if (!Invert(m_prmat)) {
341 std::cerr << m_className << "::UpdateProjectionMatrix:\n"
342 << " Inversion failed; reset to default.\n";
343 SetPlaneXY();
344 }
345}
346
347void ViewBase::Clip(const std::array<float, 3>& x0,
348 const std::array<float, 3>& x1,
349 std::array<float, 3>& xc) const {
350
351 xc.fill(0.);
352 const bool in0 = InBox(x0);
353 const bool in1 = InBox(x1);
354 if (in0 == in1) return;
355 xc = in0 ? x1 : x0;
356 const std::array<float, 3> dx = {x1[0] - x0[0], x1[1] - x0[1],
357 x1[2] - x0[2]};
358 std::array<double, 3> bmin = {m_xMinBox, m_yMinBox, m_zMinBox};
359 std::array<double, 3> bmax = {m_xMaxBox, m_yMaxBox, m_zMaxBox};
360 for (size_t i = 0; i < 3; ++i) {
361 if (dx[i] == 0. || (xc[i] >= bmin[i] && xc[i] <= bmax[i])) continue;
362 const double b = xc[i] < bmin[i] ? bmin[i] : bmax[i];
363 const double s = (b - xc[i]) / dx[i];
364 xc[i] = b;
365 for (size_t j = 0; j < 3; ++j) {
366 if (j != i) xc[j] += dx[j] * s;
367 }
368 }
369}
370
371void ViewBase::DrawLine(const std::vector<std::array<float, 3> >& xl,
372 const short col, const short lw) {
373
374 const size_t nP = xl.size();
375 if (nP < 2) return;
376
377 TGraph gr;
378 gr.SetLineColor(col);
379 gr.SetLineWidth(lw);
380
381 std::vector<float> xgr;
382 std::vector<float> ygr;
383 auto x0 = xl[0];
384 bool in0 = InBox(x0);
385 if (in0) {
386 float xp = 0., yp = 0.;
387 ToPlane(x0[0], x0[1], x0[2], xp, yp);
388 xgr.push_back(xp);
389 ygr.push_back(yp);
390 }
391 for (unsigned int j = 1; j < nP; ++j) {
392 auto x1 = xl[j];
393 bool in1 = InBox(x1);
394 if (in1 != in0) {
395 float xp = 0., yp = 0.;
396 std::array<float, 3> xc;
397 Clip(x0, x1, xc);
398 ToPlane(xc[0], xc[1], xc[2], xp, yp);
399 xgr.push_back(xp);
400 ygr.push_back(yp);
401 }
402 if (in1) {
403 float xp = 0., yp = 0.;
404 ToPlane(x1[0], x1[1], x1[2], xp, yp);
405 xgr.push_back(xp);
406 ygr.push_back(yp);
407 } else if (!xgr.empty()) {
408 gr.DrawGraph(xgr.size(), xgr.data(), ygr.data(), "Lsame");
409 xgr.clear();
410 ygr.clear();
411 }
412 x0 = x1;
413 in0 = in1;
414 }
415 if (!xgr.empty()) {
416 gr.DrawGraph(xgr.size(), xgr.data(), ygr.data(), "Lsame");
417 }
418}
419
420std::string ViewBase::LabelX() {
421
422 std::string xLabel = "";
423 constexpr double tol = 1.e-4;
424 // x portion
425 if (fabs(m_proj[0][0] - 1) < tol) {
426 xLabel = "#it{x}";
427 } else if (fabs(m_proj[0][0] + 1) < tol) {
428 xLabel = "#minus#it{x}";
429 } else if (fabs(m_proj[0][0]) > tol) {
430 xLabel = Fmt(m_proj[0][0]) + " #it{x}";
431 }
432
433 // y portion
434 if (!xLabel.empty()) {
435 if (m_proj[0][1] < -tol) {
436 xLabel += " #minus ";
437 } else if (m_proj[0][1] > tol) {
438 xLabel += " #plus ";
439 }
440 if (fabs(m_proj[0][1] - 1) < tol || fabs(m_proj[0][1] + 1) < tol) {
441 xLabel += "#it{y}";
442 } else if (fabs(m_proj[0][1]) > tol) {
443 xLabel += Fmt(fabs(m_proj[0][1])) + " #it{y}";
444 }
445 } else {
446 if (fabs(m_proj[0][1] - 1) < tol) {
447 xLabel = "#it{y}";
448 } else if (fabs(m_proj[0][1] + 1) < tol) {
449 xLabel = "#minus#it{y}";
450 } else if (fabs(m_proj[0][1]) > tol) {
451 xLabel = Fmt(m_proj[0][1]) + " #it{y}";
452 }
453 }
454
455 // z portion
456 if (!xLabel.empty()) {
457 if (m_proj[0][2] < -tol) {
458 xLabel += " #minus ";
459 } else if (m_proj[0][2] > tol) {
460 xLabel += " #plus ";
461 }
462 if (fabs(m_proj[0][2] - 1) < tol || fabs(m_proj[0][2] + 1) < tol) {
463 xLabel += "#it{z}";
464 } else if (fabs(m_proj[0][2]) > tol) {
465 xLabel += Fmt(fabs(m_proj[0][2])) + " #it{z}";
466 }
467 } else {
468 if (fabs(m_proj[0][2] - 1) < tol) {
469 xLabel = "#it{z}";
470 } else if (fabs(m_proj[0][2] + 1) < tol) {
471 xLabel = "#minus#it{z}";
472 } else if (fabs(m_proj[0][2]) > tol) {
473 xLabel = Fmt(m_proj[0][2]) + " #it{z}";
474 }
475 }
476
477 // Unit
478 xLabel += " [cm]";
479 return xLabel;
480
481}
482
483std::string ViewBase::LabelY() {
484
485 std::string yLabel = "";
486 constexpr double tol = 1.e-4;
487 // x portion
488 if (fabs(m_proj[1][0] - 1) < tol) {
489 yLabel = "#it{x}";
490 } else if (fabs(m_proj[1][0] + 1) < tol) {
491 yLabel = "#minus#it{x}";
492 } else if (fabs(m_proj[1][0]) > tol) {
493 yLabel = Fmt(m_proj[1][0]) + " #it{x}";
494 }
495
496 // y portion
497 if (!yLabel.empty()) {
498 if (m_proj[1][1] < -tol) {
499 yLabel += " #minus ";
500 } else if (m_proj[1][1] > tol) {
501 yLabel += " #plus ";
502 }
503 if (fabs(m_proj[1][1] - 1) < tol || fabs(m_proj[1][1] + 1) < tol) {
504 yLabel += "#it{y}";
505 } else if (fabs(m_proj[1][1]) > tol) {
506 yLabel += Fmt(fabs(m_proj[1][1])) + " #it{y}";
507 }
508 } else {
509 if (fabs(m_proj[1][1] - 1) < tol) {
510 yLabel = "#it{y}";
511 } else if (fabs(m_proj[1][1] + 1) < tol) {
512 yLabel = "#minus#it{y}";
513 } else if (fabs(m_proj[1][1]) > tol) {
514 yLabel = Fmt(m_proj[1][1]) + " #it{y}";
515 }
516 }
517
518 // z portion
519 if (!yLabel.empty()) {
520 if (m_proj[1][2] < -tol) {
521 yLabel += " #minus ";
522 } else if (m_proj[1][2] > tol) {
523 yLabel += " #plus ";
524 }
525 if (fabs(m_proj[1][2] - 1) < tol || fabs(m_proj[1][2] + 1) < tol) {
526 yLabel += "#it{z}";
527 } else if (fabs(m_proj[1][2]) > tol) {
528 yLabel += Fmt(fabs(m_proj[1][2])) + " #it{z}";
529 }
530 } else {
531 if (fabs(m_proj[1][2] - 1) < tol) {
532 yLabel = "#it{z}";
533 } else if (fabs(m_proj[1][2] + 1) < tol) {
534 yLabel = "#minus#it{z}";
535 } else if (fabs(m_proj[1][2]) > tol) {
536 yLabel = Fmt(m_proj[1][2]) + " #it{z}";
537 }
538 }
539
540 // Unit
541 yLabel += " [cm]";
542 return yLabel;
543}
544
545std::string ViewBase::PlaneDescription() {
546
547 std::string description;
548
549 constexpr double tol = 1.e-4;
550 // x portion
551 if (fabs(m_plane[0] - 1) < tol) {
552 description = "x";
553 } else if (fabs(m_plane[0] + 1) < tol) {
554 description = "-x";
555 } else if (fabs(m_plane[0]) > tol) {
556 description = Fmt(m_plane[0]) + " x";
557 }
558
559 // y portion
560 if (!description.empty()) {
561 if (m_plane[1] < -tol) {
562 description += " - ";
563 } else if (m_plane[1] > tol) {
564 description += " + ";
565 }
566 if (fabs(m_plane[1] - 1) < tol || fabs(m_plane[1] + 1) < tol) {
567 description += "y";
568 } else if (fabs(m_plane[1]) > tol) {
569 description += Fmt(fabs(m_plane[1])) + " y";
570 }
571 } else {
572 if (fabs(m_plane[1] - 1) < tol) {
573 description = "y";
574 } else if (fabs(m_plane[1] + 1) < tol) {
575 description = "-y";
576 } else if (fabs(m_plane[1]) > tol) {
577 description = Fmt(m_plane[1]) + " y";
578 }
579 }
580
581 // z portion
582 if (!description.empty()) {
583 if (m_plane[2] < -tol) {
584 description += " - ";
585 } else if (m_plane[2] > tol) {
586 description += " + ";
587 }
588 if (fabs(m_plane[2] - 1) < tol || fabs(m_plane[2] + 1) < tol) {
589 description += "z";
590 } else if (fabs(m_plane[2]) > tol) {
591 description += Fmt(fabs(m_plane[2])) + " z";
592 }
593 } else {
594 if (fabs(m_plane[2] - 1) < tol) {
595 description = "z";
596 } else if (fabs(m_plane[2] + 1) < tol) {
597 description = "-z";
598 } else if (fabs(m_plane[2]) > tol) {
599 description = Fmt(m_plane[2]) + " z";
600 }
601 }
602
603 // Constant
604 description += " = " + Fmt(m_plane[3]);
605 return description;
606}
607
608bool ViewBase::PlotLimits(Sensor* sensor,
609 double& xmin, double& ymin,
610 double& xmax, double& ymax) const {
611
612 if (!sensor) return false;
613 // Try to get the area/bounding box from the sensor/component.
614 std::array<double, 3> bbmin;
615 std::array<double, 3> bbmax;
616 if (!sensor->GetArea(bbmin[0], bbmin[1], bbmin[2],
617 bbmax[0], bbmax[1], bbmax[2])) {
618 std::cerr << m_className << "::PlotLimits:\n"
619 << " Sensor area is not defined.\n"
620 << " Please set the plot limits explicitly (SetArea).\n";
621 return false;
622 }
623 return PlotLimits(bbmin, bbmax, xmin, ymin, xmax, ymax);
624}
625
626bool ViewBase::PlotLimits(Component* cmp,
627 double& xmin, double& ymin,
628 double& xmax, double& ymax) const {
629
630 if (!cmp) return false;
631 // Try to get the area/bounding box from the sensor/component.
632 std::array<double, 3> bbmin;
633 std::array<double, 3> bbmax;
634 if (!cmp->GetBoundingBox(bbmin[0], bbmin[1], bbmin[2],
635 bbmax[0], bbmax[1], bbmax[2])) {
636 std::cerr << m_className << "::PlotLimits:\n"
637 << " Bounding box of the component is not defined.\n"
638 << " Please set the plot limits explicitly (SetArea).\n";
639 return false;
640 }
641 if (std::isinf(bbmin[0]) || std::isinf(bbmax[0]) ||
642 std::isinf(bbmin[1]) || std::isinf(bbmax[1]) ||
643 std::isinf(bbmin[2]) || std::isinf(bbmax[2])) {
644 std::array<double, 3> cellmin = {0., 0., 0.};
645 std::array<double, 3> cellmax = {0., 0., 0.};
646 if (!cmp->GetElementaryCell(cellmin[0], cellmin[1], cellmin[2],
647 cellmax[0], cellmax[1], cellmax[2])) {
648 std::cerr << m_className << "::PlotLimits:\n"
649 << " Cell boundaries are not defined.\n"
650 << " Please set the plot limits explicitly (SetArea).\n";
651 return false;
652 }
653 for (size_t i = 0; i < 3; ++i) {
654 if (std::isinf(bbmin[i]) || std::isinf(bbmax[i])) {
655 bbmin[i] = cellmin[i];
656 bbmax[i] = cellmax[i];
657 }
658 }
659 }
660 return PlotLimits(bbmin, bbmax, xmin, ymin, xmax, ymax);
661}
662
663bool ViewBase::PlotLimitsFromUserBox(double& xmin, double& ymin,
664 double& xmax, double& ymax) const {
665
666 std::array<double, 3> bbmin = {m_xMinBox, m_yMinBox, m_zMinBox};
667 std::array<double, 3> bbmax = {m_xMaxBox, m_yMaxBox, m_zMaxBox};
668 return PlotLimits(bbmin, bbmax, xmin, ymin, xmax, ymax);
669}
670
671bool ViewBase::PlotLimits(std::array<double, 3>& bbmin,
672 std::array<double, 3>& bbmax,
673 double& xmin, double& ymin,
674 double& xmax, double& ymax) const {
675 constexpr double tol = 1.e-4;
676 double umin[2] = {-std::numeric_limits<double>::max(),
677 -std::numeric_limits<double>::max()};
678 double umax[2] = {std::numeric_limits<double>::max(),
679 std::numeric_limits<double>::max()};
680 for (unsigned int i = 0; i < 3; ++i) {
681 bbmin[i] -= m_proj[2][i];
682 bbmax[i] -= m_proj[2][i];
683 for (unsigned int j = 0; j < 2; ++j) {
684 if (fabs(m_proj[j][i]) < tol) continue;
685 const double t1 = bbmin[i] / m_proj[j][i];
686 const double t2 = bbmax[i] / m_proj[j][i];
687 const double tmin = std::min(t1, t2);
688 const double tmax = std::max(t1, t2);
689 if (tmin > umin[j] && tmin < umax[j]) umin[j] = tmin;
690 if (tmax < umax[j] && tmax > umin[j]) umax[j] = tmax;
691 }
692 }
693 xmin = umin[0];
694 xmax = umax[0];
695 ymin = umin[1];
696 ymax = umax[1];
697 return true;
698}
699
700}
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:99
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:105
Garfield::PlottingEngine::SetDefaultStyle
void SetDefaultStyle()
Apply the default Garfield ROOT style.
Definition: PlottingEngine.cc:18
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:228
Garfield::ViewBase::m_plane
std::array< double, 4 > m_plane
Definition: ViewBase.hh:99
Garfield::ViewBase::LabelY
std::string LabelY()
Definition: ViewBase.cc:483
Garfield::ViewBase::LabelX
std::string LabelX()
Definition: ViewBase.cc:420
Garfield::ViewBase::m_prmat
std::array< std::array< double, 3 >, 3 > m_prmat
Definition: ViewBase.hh:101
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:347
Garfield::ViewBase::InBox
bool InBox(const std::array< T, 3 > &x) const
Definition: ViewBase.hh:115
Garfield::ViewBase::DrawLine
void DrawLine(const std::vector< std::array< float, 3 > > &xl, const short col, const short lw)
Definition: ViewBase.cc:371
Garfield::ViewBase::FindUnusedHistogramName
static std::string FindUnusedHistogramName(const std::string &s)
Find an unused histogram name.
Definition: ViewBase.cc:300
Garfield::ViewBase::m_className
std::string m_className
Definition: ViewBase.hh:78
Garfield::ViewBase::SetPlaneYZ
void SetPlaneYZ()
Set the viewing plane to y-z.
Definition: ViewBase.cc:284
Garfield::ViewBase::FindUnusedCanvasName
static std::string FindUnusedCanvasName(const std::string &s)
Find an unused canvas name.
Definition: ViewBase.cc:310
Garfield::ViewBase::m_proj
std::array< std::array< double, 3 >, 3 > m_proj
Definition: ViewBase.hh:96
Garfield::ViewBase::RangeSet
bool RangeSet(TPad *)
Definition: ViewBase.cc:83
Garfield::ViewBase::GetCanvas
TPad * GetCanvas()
Retrieve the canvas.
Definition: ViewBase.cc:74
Garfield::ViewBase::SetRange
void SetRange(TPad *pad, const double x0, const double y0, const double x1, const double y1)
Definition: ViewBase.cc:93
Garfield::ViewBase::SetPlaneXZ
void SetPlaneXZ()
Set the viewing plane to x-z.
Definition: ViewBase.cc:278
Garfield::ViewBase::PlotLimits
bool PlotLimits(Sensor *sensor, double &xmin, double &ymin, double &xmax, double &ymax) const
Definition: ViewBase.cc:608
Garfield::ViewBase::SetPlaneXY
void SetPlaneXY()
Set the viewing plane to x-y.
Definition: ViewBase.cc:272
Garfield::ViewBase::ToPlane
void ToPlane(const T x, const T y, const T z, T &xp, T &yp) const
Definition: ViewBase.hh:109
Garfield::ViewBase::PlaneDescription
std::string PlaneDescription()
Definition: ViewBase.cc:545
Garfield::ViewBase::PlotLimitsFromUserBox
bool PlotLimitsFromUserBox(double &xmin, double &ymin, double &xmax, double &ymax) const
Definition: ViewBase.cc:663
Garfield::ViewBase::FindUnusedFunctionName
static std::string FindUnusedFunctionName(const std::string &s)
Find an unused function name.
Definition: ViewBase.cc:290
Garfield::ViewBase::Rotate
void Rotate(const double angle)
Rotate the viewing plane (angle in radian).
Definition: ViewBase.cc:255
Garfield::ViewBase::ViewBase
ViewBase()=delete
Default constructor.
Garfield::ViewBase::UpdateProjectionMatrix
void UpdateProjectionMatrix()
Definition: ViewBase.cc:320
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