Garfield++ 3.0
A toolkit for the detailed simulation of particle detectors based on ionisation measurement in gases and semiconductors
Loading...
Searching...
No Matches
ViewIsochrons.cc
Go to the documentation of this file.
1#include <cmath>
2#include <iostream>
3#include <limits>
4#include <set>
5
6#include <TAxis.h>
7#include <TROOT.h>
8#include <TGraph.h>
9#include <TH1F.h>
10
11#include "Garfield/Sensor.hh"
14#include "Garfield/Plotting.hh"
16
17namespace {
18
19double Interpolate(const std::vector<double>& y,
20 const std::vector<double>& x, const double xx) {
21
22 const double tol = 1.e-6 * fabs(x.back() - x.front());
23 if (xx < x.front()) return y.front();
24 const auto it1 = std::upper_bound(x.cbegin(), x.cend(), xx);
25 if (it1 == x.cend()) return y.back();
26 const auto it0 = std::prev(it1);
27 const double dx = (*it1 - *it0);
28 if (dx < tol) return y[it0 - x.cbegin()];
29 const double f = (xx - *it0) / dx;
30 return y[it0 - x.cbegin()] * (1. - f) + f * y[it1 - x.cbegin()];
31}
32
33bool OnLine(const double x1, const double y1, const double x2, const double y2,
34 const double u, const double v) {
35 // Set tolerances.
36 double epsx = 1.e-10 * std::max({fabs(x1), fabs(x2), fabs(u)});
37 double epsy = 1.e-10 * std::max({fabs(y1), fabs(y2), fabs(v)});
38 epsx = std::max(1.e-10, epsx);
39 epsy = std::max(1.e-10, epsy);
40
41 if ((fabs(x1 - u) <= epsx && fabs(y1 - v) <= epsy) ||
42 (fabs(x2 - u) <= epsx && fabs(y2 - v) <= epsy)) {
43 // Point to be examined coincides with start or end.
44 return true;
45 } else if (fabs(x1 - x2) <= epsx && fabs(y1 - y2) <= epsy) {
46 // The line (x1, y1) to (x2, y2) is in fact a point.
47 return false;
48 }
49 double xc = 0., yc = 0.;
50 if (fabs(u - x1) + fabs(v - y1) < fabs(u - x2) + fabs(v - y2)) {
51 // (u, v) is nearer to (x1, y1).
52 const double dx = (x2 - x1);
53 const double dy = (y2 - y1);
54 const double xl = ((u - x1) * dx + (v - y1) * dy) / (dx * dx + dy * dy);
55 if (xl < 0.) {
56 xc = x1;
57 yc = y1;
58 } else if (xl > 1.) {
59 xc = x2;
60 yc = y2;
61 } else {
62 xc = x1 + xl * dx;
63 yc = y1 + xl * dy;
64 }
65 } else {
66 // (u, v) is nearer to (x2, y2).
67 const double dx = (x1 - x2);
68 const double dy = (y1 - y2);
69 const double xl = ((u - x2) * dx + (v - y2) * dy) / (dx * dx + dy * dy);
70 if (xl < 0.) {
71 xc = x2;
72 yc = y2;
73 } else if (xl > 1.) {
74 xc = x1;
75 yc = y1;
76 } else {
77 xc = x2 + xl * dx;
78 yc = y2 + xl * dy;
79 }
80 }
81 // See whether the point is on the line.
82 if (fabs(u - xc) < epsx && fabs(v - yc) < epsy) {
83 return true;
84 }
85 return false;
86}
87
88bool Crossing(const double x1, const double y1, const double x2,
89 const double y2, const double u1, const double v1,
90 const double u2, const double v2) {
91
92 // Check for a point of one line located on the other line.
93 if (OnLine(x1, y1, x2, y2, u1, v1) || OnLine(x1, y1, x2, y2, u2, v2) ||
94 OnLine(u1, v1, u2, v2, x1, y1) || OnLine(u1, v1, u2, v2, x2, y2)) {
95 return true;
96 }
97 // Matrix to compute the crossing point.
98 std::array<std::array<double, 2>, 2> a;
99 a[0][0] = y2 - y1;
100 a[0][1] = v2 - v1;
101 a[1][0] = x1 - x2;
102 a[1][1] = u1 - u2;
103 const double det = a[0][0] * a[1][1] - a[1][0] * a[0][1];
104 // Set tolerances.
105 double epsx = 1.e-10 * std::max({fabs(x1), fabs(x2), fabs(u1), fabs(u2)});
106 double epsy = 1.e-10 * std::max({fabs(y1), fabs(y2), fabs(v1), fabs(v2)});
107 epsx = std::max(epsx, 1.e-10);
108 epsy = std::max(epsy, 1.e-10);
109 if (fabs(det) < epsx * epsy) {
110 // Parallel, non-touching.
111 return false;
112 }
113 // Crossing, non-trivial lines: solve crossing equations.
114 const double aux = a[1][1];
115 a[1][1] = a[0][0] / det;
116 a[0][0] = aux / det;
117 a[1][0] = -a[1][0] / det;
118 a[0][1] = -a[0][1] / det;
119 // Compute crossing point.
120 const double xc = a[0][0] * (x1 * y2 - x2 * y1) + a[1][0] * (u1 * v2 - u2 * v1);
121 const double yc = a[0][1] * (x1 * y2 - x2 * y1) + a[1][1] * (u1 * v2 - u2 * v1);
122 // See whether the crossing point is on both lines.
123 if (OnLine(x1, y1, x2, y2, xc, yc) && OnLine(u1, v1, u2, v2, xc, yc)) {
124 // Intersecting lines.
125 return true;
126 }
127 // Crossing point not on both lines.
128 return false;
129}
130
131}
132
133namespace Garfield {
134
137}
138
140 if (!s) {
141 std::cerr << m_className << "::SetSensor: Null pointer.\n";
142 return;
143 }
144
145 m_sensor = s;
146 m_component = nullptr;
147}
148
150 if (!c) {
151 std::cerr << m_className << "::SetComponent: Null pointer.\n";
152 return;
153 }
154
155 m_component = c;
156 m_sensor = nullptr;
157}
158
159void ViewIsochrons::SetArea(const double xmin, const double ymin,
160 const double xmax, const double ymax) {
161 // Check range, assign if non-null.
162 if (xmin == xmax || ymin == ymax) {
163 std::cerr << m_className << "::SetArea: Null area range is not permitted.\n"
164 << " " << xmin << " < x < " << xmax << "\n"
165 << " " << ymin << " < y < " << ymax << "\n";
166 return;
167 }
168 m_xmin = std::min(xmin, xmax);
169 m_ymin = std::min(ymin, ymax);
170 m_xmax = std::max(xmin, xmax);
171 m_ymax = std::max(ymin, ymax);
172 m_hasUserArea = true;
173}
174
176
177 if (ar < 0.) {
178 std::cerr << m_className << "::SetAspectRatioSwitch: Value must be > 0.\n";
179 return;
180 }
181 m_aspectRatio = ar;
182}
183
184void ViewIsochrons::SetLoopThreshold(const double thr) {
185
186 if (thr < 0. || thr > 1.) {
187 std::cerr << m_className << "::SetLoopThreshold:\n"
188 << " Value must be between 0 and 1.\n";
189 return;
190 }
191 m_loopThreshold = thr;
192}
193
195
196 if (thr < 0. || thr > 1.) {
197 std::cerr << m_className << "::SetConnectionThreshold:\n"
198 << " Value must be between 0 and 1.\n";
199 return;
200 }
201 m_connectionThreshold = thr;
202}
203
204void ViewIsochrons::PlotIsochrons(const double tstep,
205 const std::vector<std::array<double, 3> >& points, const bool rev,
206 const bool colour, const bool markers, const bool plotDriftLines) {
207
208 if (!m_sensor && !m_component) {
209 std::cerr << m_className << "::PlotIsochrons:\n"
210 << " Neither sensor nor component are defined.\n";
211 return;
212 }
213 if (tstep <= 0.) {
214 std::cerr << m_className << "::PlotIsochrons: Time step must be > 0.\n";
215 return;
216 }
217 if (points.empty()) {
218 std::cerr << m_className << "::PlotIsochrons:\n"
219 << " No starting points provided.\n";
220 return;
221 }
222 SetupCanvas();
223 if (!Range()) return;
224 auto frame = m_canvas->DrawFrame(m_xmin, m_ymin, m_xmax, m_ymax);
225 frame->GetXaxis()->SetTitle(m_xLabel);
226 frame->GetYaxis()->SetTitle(m_yLabel);
227 m_canvas->Update();
228
229 //-----------------------------------------------------------------------
230 // DRFEQT - The main routine (DRFEQT) accumulates equal drift time data
231 // DRFEQP which is plotted as a set of contours in the entry DRFEQP.
232 //-----------------------------------------------------------------------
233 std::vector<std::vector<std::array<double, 3> > > driftLines;
234 std::vector<std::array<double, 3> > startPoints;
235 std::vector<std::array<double, 3> > endPoints;
236 std::vector<int> statusCodes;
237 // Accumulate drift lines.
238 ComputeDriftLines(tstep, points, driftLines, startPoints, endPoints,
239 statusCodes, rev);
240 const unsigned int nDriftLines = driftLines.size();
241 if (nDriftLines < 2) {
242 std::cerr << m_className << "::PlotIsochrons: Too few drift lines.\n";
243 return;
244 }
245 // Keep track of the largest number of contours.
246 std::size_t nContours = 0;
247 for (const auto& driftLine : driftLines) {
248 nContours = std::max(nContours, driftLine.size());
249 }
250 if (nContours == 0) {
251 std::cerr << m_className << "::PlotIsochrons: No contour lines.\n";
252 return;
253 }
254
255 std::set<int> allStats;
256 for (const auto stat : statusCodes) allStats.insert(stat);
257
258 // DRFEQP
259 if (m_debug) {
260 std::cout << m_className << "::PlotIsochrons:\n"
261 << " Drawing " << nContours << " contours, "
262 << nDriftLines << " drift lines.\n";
263 std::printf(" Connection threshold: %10.3f\n", m_connectionThreshold);
264 std::printf(" Aspect ratio threshold: %10.3f\n", m_aspectRatio);
265 std::printf(" Loop closing threshold: %10.3f\n", m_loopThreshold);
266 if (m_sortContours) {
267 std::cout << " Sort contours.\n";
268 } else {
269 std::cout << " Do not sort contours.\n";
270 }
271 if (m_checkCrossings) {
272 std::cout << " Check for crossings.\n";
273 } else {
274 std::cout << " Do not check for crossings.\n";
275 }
276 if (markers) {
277 std::cout << " Mark isochron points.\n";
278 } else {
279 std::cout << " Draw isochron lines.\n";
280 }
281 }
282
283 // Loop over the equal time contours.
284 TGraph graph;
285 graph.SetLineColor(kGray + 2);
286 graph.SetMarkerColor(kGray + 2);
287 graph.SetLineStyle(m_lineStyle);
288 graph.SetMarkerStyle(m_markerStyle);
289 const double colRange = double(gStyle->GetNumberOfColors()) / nContours;
290 for (unsigned int ic = 0; ic < nContours; ++ic) {
291 if (colour) {
292 const auto col = gStyle->GetColorPalette(int((ic + 0.99) * colRange));
293 graph.SetLineColor(col);
294 graph.SetMarkerColor(col);
295 }
296 for (int stat : allStats) {
297 std::vector<std::pair<std::array<double, 4>, unsigned int> > contour;
298 // Loop over the drift lines, picking up the points when OK.
299 for (unsigned int k = 0; k < nDriftLines; ++k) {
300 const auto& dl = driftLines[k];
301 // Reject any undesirable combinations.
302 if (statusCodes[k] != stat || ic >= dl.size()) continue;
303 // Add the point to the contour line.
304 std::array<double, 4> point = {dl[ic][0], dl[ic][1], dl[ic][2], 0.};
305 contour.push_back(std::make_pair(point, k));
306 }
307 // Skip the plot of this contour if there are no points.
308 if (contour.empty()) continue;
309 bool circle = false;
310 // If requested, sort the points on the contour line.
311 if (m_sortContours && !markers) SortContour(contour, circle);
312 // Plot this contour.
313 if (markers) {
314 // Simply mark the contours if this was requested.
315 std::vector<double> xp;
316 std::vector<double> yp;
317 std::vector<double> zp;
318 for (const auto& point : contour) {
319 const double x = point.first[0];
320 const double y = point.first[1];
321 const double z = point.first[2];
322 xp.push_back(m_proj[0][0] * x + m_proj[1][0] * y + z * m_plane[0]);
323 yp.push_back(m_proj[0][1] * x + m_proj[1][1] * y + z * m_plane[1]);
324 zp.push_back(m_proj[0][2] * x + m_proj[1][2] * y + z * m_plane[2]);
325 }
326 graph.DrawGraph(xp.size(), xp.data(), yp.data(), "Psame");
327 continue;
328 }
329 // Regular plotting.
330 const double tolx = (m_xmax - m_xmin) * m_connectionThreshold;
331 const double toly = (m_ymax - m_ymin) * m_connectionThreshold;
332 // Flag to keep track if the segment is interrupted by a drift line
333 // or if it is too long.
334 bool gap = false;
335 // Coordinates to be plotted.
336 std::vector<double> xp;
337 std::vector<double> yp;
338 std::vector<double> zp;
339 const unsigned int nP = contour.size();
340 for (unsigned int i = 0; i < nP; ++i) {
341 gap = false;
342 const auto x0 = contour[i].first[0];
343 const auto y0 = contour[i].first[1];
344 const auto z0 = contour[i].first[2];
345 xp.push_back(m_proj[0][0] * x0 + m_proj[1][0] * y0 + z0 * m_plane[0]);
346 yp.push_back(m_proj[0][1] * x0 + m_proj[1][1] * y0 + z0 * m_plane[1]);
347 zp.push_back(m_proj[0][2] * x0 + m_proj[1][2] * y0 + z0 * m_plane[2]);
348 if (i == nP - 1) break;
349 const auto x1 = contour[i + 1].first[0];
350 const auto y1 = contour[i + 1].first[1];
351 // Reject contour segments which are long compared with AREA.
352 if (fabs(x1 - x0) > tolx || fabs(y1 - y0) > toly) gap = true;
353 // Get the indices of the drift lines corresponding
354 // to these two points on the contour line.
355 const auto i0 = contour[i].second;
356 const auto i1 = contour[i + 1].second;
357 // Set the BREAK flag if it crosses some stored drift line segment.
358 if (m_checkCrossings && !gap) {
359 for (unsigned int k = 0; k < nDriftLines; ++k) {
360 const auto& dl = driftLines[k];
361 for (unsigned int jc = 0; jc < dl.size(); ++jc) {
362 if ((i0 == k || i1 == k) && (jc == ic || jc + 1 == ic)) {
363 continue;
364 }
365 const auto& p0 = dl[jc];
366 const auto& p1 = jc == dl.size() - 1 ? endPoints[k] : dl[jc + 1];
367 if (Crossing(p0[0], p0[1], p1[0], p1[1], x0, y0, x1, y1)) {
368 gap = true;
369 break;
370 }
371 }
372 if (gap) break;
373 if ((i0 == k || i1 == k) && ic == 0) continue;
374 const auto& p0 = startPoints[k];
375 if (Crossing(p0[0], p0[1], dl[0][0], dl[0][1],
376 x0, y0, x1, y1)) {
377 gap = true;
378 break;
379 }
380 }
381 }
382 // If there has been a break, plot what we have already.
383 if (gap) {
384 if (xp.size() > 1) {
385 // Plot line.
386 graph.DrawGraph(xp.size(), xp.data(), yp.data(), "Lsame");
387 } else {
388 // Plot markers.
389 graph.DrawGraph(xp.size(), xp.data(), yp.data(), "Psame");
390 }
391 xp.clear();
392 yp.clear();
393 zp.clear();
394 }
395 }
396 // Plot the remainder.
397 if (xp.size() > 1) {
398 graph.DrawGraph(xp.size(), xp.data(), yp.data(), "Lsame");
399 } else if (!xp.empty()) {
400 graph.DrawGraph(xp.size(), xp.data(), yp.data(), "Psame");
401 }
402 }
403 }
404
405 gPad->Update();
406 if (!plotDriftLines) return;
407
408 graph.SetLineStyle(1);
409 if (m_particle == Particle::Electron) {
410 graph.SetLineColor(kOrange - 3);
411 } else {
412 graph.SetLineColor(kRed + 1);
413 }
414 for (unsigned int i = 0; i < nDriftLines; ++i) {
415 std::vector<double> xp;
416 std::vector<double> yp;
417 const double x0 = startPoints[i][0];
418 const double y0 = startPoints[i][1];
419 const double z0 = startPoints[i][2];
420 xp.push_back(m_proj[0][0] * x0 + m_proj[1][0] * y0 + z0 * m_plane[0]);
421 yp.push_back(m_proj[0][1] * x0 + m_proj[1][1] * y0 + z0 * m_plane[1]);
422 for (const auto& point : driftLines[i]) {
423 const double x = point[0];
424 const double y = point[1];
425 const double z = point[2];
426 xp.push_back(m_proj[0][0] * x + m_proj[1][0] * y + z * m_plane[0]);
427 yp.push_back(m_proj[0][1] * x + m_proj[1][1] * y + z * m_plane[1]);
428 }
429 const double x1 = endPoints[i][0];
430 const double y1 = endPoints[i][1];
431 const double z1 = endPoints[i][2];
432 xp.push_back(m_proj[0][0] * x1 + m_proj[1][0] * y1 + z1 * m_plane[0]);
433 yp.push_back(m_proj[0][1] * x1 + m_proj[1][1] * y1 + z1 * m_plane[1]);
434 graph.DrawGraph(xp.size(), xp.data(), yp.data(), "Lsame");
435 }
436}
437
438void ViewIsochrons::ComputeDriftLines(const double tstep,
439 const std::vector<std::array<double, 3> >& points,
440 std::vector<std::vector<std::array<double, 3> > >& driftLines,
441 std::vector<std::array<double, 3> >& startPoints,
442 std::vector<std::array<double, 3> >& endPoints,
443 std::vector<int>& statusCodes, const bool rev) {
444
445 DriftLineRKF drift;
446 Sensor sensor;
447 if (m_sensor) {
448 drift.SetSensor(m_sensor);
449 } else {
450 sensor.AddComponent(m_component);
451 drift.SetSensor(&sensor);
452 }
453 const double lx = 0.1 * fabs(m_xmax - m_xmin);
454 const double ly = 0.1 * fabs(m_ymax - m_ymin);
455 drift.SetMaximumStepSize(std::min(lx, ly));
456 drift.EnableSignalCalculation(false);
457 for (const auto& point : points) {
458 if (m_particle == Particle::Electron) {
459 if (m_positive) {
460 drift.DriftPositron(point[0], point[1], point[2], 0.);
461 } else {
462 drift.DriftElectron(point[0], point[1], point[2], 0.);
463 }
464 } else {
465 if (m_positive) {
466 drift.DriftIon(point[0], point[1], point[2], 0.);
467 } else {
468 drift.DriftNegativeIon(point[0], point[1], point[2], 0.);
469 }
470 }
471 const unsigned int nu = drift.GetNumberOfDriftLinePoints();
472 // Check that the drift line has enough points.
473 if (nu < 3) continue;
474 int status = 0;
475 double xf = 0., yf = 0., zf = 0., tf = 0.;
476 drift.GetEndPoint(xf, yf, zf, tf, status);
477 // Find the number of points to be stored.
478 const unsigned int nSteps = static_cast<unsigned int>(tf / tstep);
479 if (nSteps == 0) continue;
480 std::vector<double> xu(nu, 0.);
481 std::vector<double> yu(nu, 0.);
482 std::vector<double> zu(nu, 0.);
483 std::vector<double> tu(nu, 0.);
484 for (unsigned int i = 0; i < nu; ++i) {
485 drift.GetDriftLinePoint(i, xu[i], yu[i], zu[i], tu[i]);
486 }
487 if (rev) {
488 for (auto& t : tu) t = tf - t;
489 std::reverse(std::begin(xu), std::end(xu));
490 std::reverse(std::begin(yu), std::end(yu));
491 std::reverse(std::begin(zu), std::end(zu));
492 std::reverse(std::begin(tu), std::end(tu));
493 }
494 std::vector<std::array<double, 3> > tab;
495 // Interpolate at regular time intervals.
496 for (unsigned int i = 0; i < nSteps; ++i) {
497 const double t = (i + 1) * tstep;
498 // tab.push_back(PLACO3(Interpolate(xu, tu, t),
499 // Interpolate(yu, tu, t),
500 // Interpolate(zu, tu, t)));
501 std::array<double, 3> step = {Interpolate(xu, tu, t),
502 Interpolate(yu, tu, t),
503 Interpolate(zu, tu, t)};
504 tab.push_back(step);
505 }
506 driftLines.push_back(std::move(tab));
507 std::array<double, 3> start = {xu[0], yu[0], zu[0]};
508 std::array<double, 3> end = {xu[nu - 1], yu[nu - 1], zu[nu - 1]};
509 startPoints.push_back(std::move(start));
510 endPoints.push_back(std::move(end));
511 // Store the drift line return code.
512 if (rev) {
513 statusCodes.push_back(status);
514 } else {
515 statusCodes.push_back(0);
516 }
517 }
518}
519
521 // Default projection: x-y at z=0
522 m_proj[0][0] = 1;
523 m_proj[1][0] = 0;
524 m_proj[2][0] = 0;
525 m_proj[0][1] = 0;
526 m_proj[1][1] = 1;
527 m_proj[2][1] = 0;
528 m_proj[0][2] = 0;
529 m_proj[1][2] = 0;
530 m_proj[2][2] = 0;
531
532 // Plane description
533 m_plane[0] = 0;
534 m_plane[1] = 0;
535 m_plane[2] = 1;
536 m_plane[3] = 0;
537
538 // Prepare axis labels.
539 Labels();
540}
541
542void ViewIsochrons::Labels() {
543 // Initialisation of the x-axis label
544 strcpy(m_xLabel, "\0");
545 char buf[100];
546
547 const double tol = 1.e-4;
548 // x portion
549 if (fabs(m_proj[0][0] - 1) < tol) {
550 strcat(m_xLabel, "x");
551 } else if (fabs(m_proj[0][0] + 1) < tol) {
552 strcat(m_xLabel, "-x");
553 } else if (m_proj[0][0] > tol) {
554 sprintf(buf, "%g x", m_proj[0][0]);
555 strcat(m_xLabel, buf);
556 } else if (m_proj[0][0] < -tol) {
557 sprintf(buf, "%g x", m_proj[0][0]);
558 strcat(m_xLabel, buf);
559 }
560
561 // y portion
562 if (strlen(m_xLabel) > 0) {
563 if (m_proj[0][1] < -tol) {
564 strcat(m_xLabel, " - ");
565 } else if (m_proj[0][1] > tol) {
566 strcat(m_xLabel, " + ");
567 }
568 if (fabs(m_proj[0][1] - 1) < tol || fabs(m_proj[0][1] + 1) < tol) {
569 strcat(m_xLabel, "y");
570 } else if (fabs(m_proj[0][1]) > tol) {
571 sprintf(buf, "%g y", fabs(m_proj[0][1]));
572 strcat(m_xLabel, buf);
573 }
574 } else {
575 if (fabs(m_proj[0][1] - 1) < tol) {
576 strcat(m_xLabel, "y");
577 } else if (fabs(m_proj[0][1] + 1) < tol) {
578 strcat(m_xLabel, "-y");
579 } else if (m_proj[0][1] > tol) {
580 sprintf(buf, "%g y", m_proj[0][1]);
581 strcat(m_xLabel, buf);
582 } else if (m_proj[0][1] < -tol) {
583 sprintf(buf, "%g y", m_proj[0][1]);
584 strcat(m_xLabel, buf);
585 }
586 }
587
588 // z portion
589 if (strlen(m_xLabel) > 0) {
590 if (m_proj[0][2] < -tol) {
591 strcat(m_xLabel, " - ");
592 } else if (m_proj[0][2] > tol) {
593 strcat(m_xLabel, " + ");
594 }
595 if (fabs(m_proj[0][2] - 1) < tol || fabs(m_proj[0][2] + 1) < tol) {
596 strcat(m_xLabel, "z");
597 } else if (fabs(m_proj[0][2]) > tol) {
598 sprintf(buf, "%g z", fabs(m_proj[0][2]));
599 strcat(m_xLabel, buf);
600 }
601 } else {
602 if (fabs(m_proj[0][2] - 1) < tol) {
603 strcat(m_xLabel, "z");
604 } else if (fabs(m_proj[0][2] + 1) < tol) {
605 strcat(m_xLabel, "-z");
606 } else if (m_proj[0][2] > tol) {
607 sprintf(buf, "%g z", m_proj[0][2]);
608 strcat(m_xLabel, buf);
609 } else if (m_proj[0][2] < -tol) {
610 sprintf(buf, "%g z", m_proj[0][2]);
611 strcat(m_xLabel, buf);
612 }
613 }
614
615 // Unit
616 strcat(m_xLabel, " [cm]");
617
618 // Initialisation of the y-axis label
619 strcpy(m_yLabel, "\0");
620
621 // x portion
622 if (fabs(m_proj[1][0] - 1) < tol) {
623 strcat(m_yLabel, "x");
624 } else if (fabs(m_proj[1][0] + 1) < tol) {
625 strcat(m_yLabel, "-x");
626 } else if (m_proj[1][0] > tol) {
627 sprintf(buf, "%g x", m_proj[1][0]);
628 strcat(m_yLabel, buf);
629 } else if (m_proj[1][0] < -tol) {
630 sprintf(buf, "%g x", m_proj[1][0]);
631 strcat(m_yLabel, buf);
632 }
633
634 // y portion
635 if (strlen(m_yLabel) > 0) {
636 if (m_proj[1][1] < -tol) {
637 strcat(m_yLabel, " - ");
638 } else if (m_proj[1][1] > tol) {
639 strcat(m_yLabel, " + ");
640 }
641 if (fabs(m_proj[1][1] - 1) < tol || fabs(m_proj[1][1] + 1) < tol) {
642 strcat(m_yLabel, "y");
643 } else if (fabs(m_proj[1][1]) > tol) {
644 sprintf(buf, "%g y", fabs(m_proj[1][1]));
645 strcat(m_yLabel, buf);
646 }
647 } else {
648 if (fabs(m_proj[1][1] - 1) < tol) {
649 strcat(m_yLabel, "y");
650 } else if (fabs(m_proj[1][1] + 1) < tol) {
651 strcat(m_yLabel, "-y");
652 } else if (m_proj[1][1] > tol) {
653 sprintf(buf, "%g y", m_proj[1][1]);
654 strcat(m_yLabel, buf);
655 } else if (m_proj[1][1] < -tol) {
656 sprintf(buf, "%g y", m_proj[1][1]);
657 strcat(m_yLabel, buf);
658 }
659 }
660
661 // z portion
662 if (strlen(m_yLabel) > 0) {
663 if (m_proj[1][2] < -tol) {
664 strcat(m_yLabel, " - ");
665 } else if (m_proj[1][2] > tol) {
666 strcat(m_yLabel, " + ");
667 }
668 if (fabs(m_proj[1][2] - 1) < tol || fabs(m_proj[1][2] + 1) < tol) {
669 strcat(m_yLabel, "z");
670 } else if (fabs(m_proj[1][2]) > tol) {
671 sprintf(buf, "%g z", fabs(m_proj[1][2]));
672 strcat(m_yLabel, buf);
673 }
674 } else {
675 if (fabs(m_proj[1][2] - 1) < tol) {
676 strcat(m_yLabel, "z");
677 } else if (fabs(m_proj[1][2] + 1) < tol) {
678 strcat(m_yLabel, "-z");
679 } else if (m_proj[1][2] > tol) {
680 sprintf(buf, "%g z", m_proj[1][2]);
681 strcat(m_yLabel, buf);
682 } else if (m_proj[1][2] < -tol) {
683 sprintf(buf, "%g z", m_proj[1][2]);
684 strcat(m_yLabel, buf);
685 }
686 }
687
688 // Unit
689 strcat(m_yLabel, " [cm]");
690
691 // Initialisation of the plane label
692 strcpy(m_description, "\0");
693
694 // x portion
695 if (fabs(m_plane[0] - 1) < tol) {
696 strcat(m_description, "x");
697 } else if (fabs(m_plane[0] + 1) < tol) {
698 strcat(m_description, "-x");
699 } else if (m_plane[0] > tol) {
700 sprintf(buf, "%g x", m_plane[0]);
701 strcat(m_description, buf);
702 } else if (m_plane[0] < -tol) {
703 sprintf(buf, "%g x", m_plane[0]);
704 strcat(m_description, buf);
705 }
706
707 // y portion
708 if (strlen(m_description) > 0) {
709 if (m_plane[1] < -tol) {
710 strcat(m_description, " - ");
711 } else if (m_plane[1] > tol) {
712 strcat(m_description, " + ");
713 }
714 if (fabs(m_plane[1] - 1) < tol || fabs(m_plane[1] + 1) < tol) {
715 strcat(m_description, "y");
716 } else if (fabs(m_plane[1]) > tol) {
717 sprintf(buf, "%g y", fabs(m_plane[1]));
718 strcat(m_description, buf);
719 }
720 } else {
721 if (fabs(m_plane[1] - 1) < tol) {
722 strcat(m_description, "y");
723 } else if (fabs(m_plane[1] + 1) < tol) {
724 strcat(m_description, "-y");
725 } else if (m_plane[1] > tol) {
726 sprintf(buf, "%g y", m_plane[1]);
727 strcat(m_description, buf);
728 } else if (m_plane[1] < -tol) {
729 sprintf(buf, "%g y", m_plane[1]);
730 strcat(m_description, buf);
731 }
732 }
733
734 // z portion
735 if (strlen(m_description) > 0) {
736 if (m_plane[2] < -tol) {
737 strcat(m_description, " - ");
738 } else if (m_plane[2] > tol) {
739 strcat(m_description, " + ");
740 }
741 if (fabs(m_plane[2] - 1) < tol || fabs(m_plane[2] + 1) < tol) {
742 strcat(m_description, "z");
743 } else if (fabs(m_plane[2]) > tol) {
744 sprintf(buf, "%g z", fabs(m_plane[2]));
745 strcat(m_description, buf);
746 }
747 } else {
748 if (fabs(m_plane[2] - 1) < tol) {
749 strcat(m_description, "z");
750 } else if (fabs(m_plane[2] + 1) < tol) {
751 strcat(m_description, "-z");
752 } else if (m_plane[2] > tol) {
753 sprintf(buf, "%g z", m_plane[2]);
754 strcat(m_description, buf);
755 } else if (m_plane[2] < -tol) {
756 sprintf(buf, "%g z", m_plane[2]);
757 strcat(m_description, buf);
758 }
759 }
760
761 // Constant
762 sprintf(buf, " = %g", m_plane[3]);
763 strcat(m_description, buf);
764
765 if (m_debug) {
766 std::cout << m_className << "::Labels:\n"
767 << " x label: |" << m_xLabel << "|\n"
768 << " y label: |" << m_yLabel << "|\n"
769 << " plane: |" << m_description << "|\n";
770 }
771}
772
773void ViewIsochrons::SetPlane(const double fx, const double fy, const double fz,
774 const double x0, const double y0, const double z0) {
775 // Calculate two in-plane vectors for the normal vector
776 const double fnorm = sqrt(fx * fx + fy * fy + fz * fz);
777 if (fnorm > 0 && fx * fx + fz * fz > 0) {
778 const double fxz = sqrt(fx * fx + fz * fz);
779 m_proj[0][0] = fz / fxz;
780 m_proj[0][1] = 0;
781 m_proj[0][2] = -fx / fxz;
782 m_proj[1][0] = -fx * fy / (fxz * fnorm);
783 m_proj[1][1] = (fx * fx + fz * fz) / (fxz * fnorm);
784 m_proj[1][2] = -fy * fz / (fxz * fnorm);
785 m_proj[2][0] = x0;
786 m_proj[2][1] = y0;
787 m_proj[2][2] = z0;
788 } else if (fnorm > 0 && fy * fy + fz * fz > 0) {
789 const double fyz = sqrt(fy * fy + fz * fz);
790 m_proj[0][0] = (fy * fy + fz * fz) / (fyz * fnorm);
791 m_proj[0][1] = -fx * fz / (fyz * fnorm);
792 m_proj[0][2] = -fy * fz / (fyz * fnorm);
793 m_proj[1][0] = 0;
794 m_proj[1][1] = fz / fyz;
795 m_proj[1][2] = -fy / fyz;
796 m_proj[2][0] = x0;
797 m_proj[2][1] = y0;
798 m_proj[2][2] = z0;
799 } else {
800 std::cout << m_className << "::SetPlane:\n"
801 << " Normal vector has zero norm. No new projection set.\n";
802 }
803
804 // Store the plane description
805 m_plane[0] = fx;
806 m_plane[1] = fy;
807 m_plane[2] = fz;
808 m_plane[3] = fx * x0 + fy * y0 + fz * z0;
809
810 // Make labels to be placed along the axes
811 Labels();
812}
813
814void ViewIsochrons::Rotate(const double theta) {
815 // Rotate the axes
816 double auxu[3], auxv[3];
817 const double ctheta = cos(theta);
818 const double stheta = sin(theta);
819 for (int i = 0; i < 3; ++i) {
820 auxu[i] = ctheta * m_proj[0][i] - stheta * m_proj[1][i];
821 auxv[i] = stheta * m_proj[0][i] + ctheta * m_proj[1][i];
822 }
823 for (int i = 0; i < 3; ++i) {
824 m_proj[0][i] = auxu[i];
825 m_proj[1][i] = auxv[i];
826 }
827
828 // Make labels to be placed along the axes
829 Labels();
830}
831
832void ViewIsochrons::SetupCanvas() {
833 if (!m_canvas) {
834 m_canvas = new TCanvas();
835 m_canvas->SetTitle("Isochrons");
836 m_hasExternalCanvas = false;
837 }
838 m_canvas->cd();
839}
840
841bool ViewIsochrons::Range() {
842
843 if (m_hasUserArea) return true;
844 // Try to get the area/bounding box from the sensor/component.
845 double bbmin[3];
846 double bbmax[3];
847 if (m_sensor) {
848 if (!m_sensor->GetArea(bbmin[0], bbmin[1], bbmin[2], bbmax[0], bbmax[1],
849 bbmax[2])) {
850 std::cerr << m_className << "::Range:\n"
851 << " Sensor area is not defined.\n"
852 << " Please set the plot range explicitly (SetArea).\n";
853 return false;
854 }
855 } else {
856 if (!m_component->GetBoundingBox(bbmin[0], bbmin[1], bbmin[2], bbmax[0],
857 bbmax[1], bbmax[2])) {
858 std::cerr << m_className << "::Range:\n"
859 << " Bounding box of the component is not defined.\n"
860 << " Please set the plot range explicitly (SetArea).\n";
861 return false;
862 }
863 }
864 const double tol = 1.e-4;
865 double umin[2] = {-std::numeric_limits<double>::max(),
866 -std::numeric_limits<double>::max()};
867 double umax[2] = {std::numeric_limits<double>::max(),
868 std::numeric_limits<double>::max()};
869 for (unsigned int i = 0; i < 3; ++i) {
870 bbmin[i] -= m_proj[2][i];
871 bbmax[i] -= m_proj[2][i];
872 for (unsigned int j = 0; j < 2; ++j) {
873 if (fabs(m_proj[j][i]) < tol) continue;
874 const double t1 = bbmin[i] / m_proj[j][i];
875 const double t2 = bbmax[i] / m_proj[j][i];
876 const double tmin = std::min(t1, t2);
877 const double tmax = std::max(t1, t2);
878 if (tmin > umin[j] && tmin < umax[j]) umin[j] = tmin;
879 if (tmax < umax[j] && tmax > umin[j]) umax[j] = tmax;
880 }
881 }
882 m_xmin = umin[0];
883 m_xmax = umax[0];
884 m_ymin = umin[1];
885 m_ymax = umax[1];
886 return true;
887}
888
889void ViewIsochrons::SortContour(
890 std::vector<std::pair<std::array<double, 4>, unsigned int> >& contour,
891 bool& circle) {
892
893 if (contour.size() < 2) return;
894 // First compute the centre of gravity.
895 double xcog = 0.;
896 double ycog = 0.;
897 for (const auto& point : contour) {
898 xcog += point.first[0];
899 ycog += point.first[1];
900 }
901 const double scale = 1. / contour.size();
902 xcog *= scale;
903 ycog *= scale;
904 // Compute angles wrt to the centre of gravity and principal axes.
905 double sxx = 0.;
906 double sxy = 0.;
907 double syy = 0.;
908 for (const auto& point : contour) {
909 const double dx = point.first[0] - xcog;
910 const double dy = point.first[1] - ycog;
911 sxx += dx * dx;
912 sxy += dx * dy;
913 syy += dy * dy;
914 }
915 const double theta = 0.5 * atan2(2 * sxy, sxx - syy);
916 const double ct = cos(theta);
917 const double st = sin(theta);
918 // Evaluate dispersions around the principal axes.
919 sxx = 0.;
920 syy = 0.;
921 for (const auto& point : contour) {
922 const double dx = point.first[0] - xcog;
923 const double dy = point.first[1] - ycog;
924 sxx += fabs(+ct * dx + st * dy);
925 syy += fabs(-st * dx + ct * dy);
926 }
927 // Decide whether this is more linear or more circular.
928 if (fabs(sxx) > m_aspectRatio * fabs(syy) ||
929 fabs(syy) > m_aspectRatio * fabs(sxx)) {
930 circle = false;
931 } else {
932 circle = true;
933 }
934 // Set a sorting coordinate accordingly.
935 for (auto& point : contour) {
936 const double dx = point.first[0] - xcog;
937 const double dy = point.first[1] - ycog;
938 point.first[3] = circle ? atan2(dy, dx) : ct * dx + st * dy;
939 }
940 // Sort the points.
941 std::sort(contour.begin(), contour.end(),
942 [](const std::pair<std::array<double, 4>, int>& p1,
943 const std::pair<std::array<double, 4>, int>& p2) {
944 return p1.first[3] < p2.first[3]; }
945 );
946 if (!circle) return;
947 // For circles, perhaps add the first point to the end of the list.
948 // Compute breakpoint, total distance and maximum distance.
949 double dsum = 0.;
950 double dmax = -1.;
951 unsigned int imax = 0;
952 const unsigned int nPoints = contour.size();
953 for (unsigned int j = 0; j < nPoints; ++j) {
954 const auto& p1 = contour[j];
955 const auto& p0 = j > 0 ? contour[j - 1] : contour.back();
956 const double dx = p1.first[0] - p0.first[0];
957 const double dy = p1.first[1] - p0.first[1];
958 const double d = sqrt(dx * dx + dy * dy);
959 if (j > 0) dsum += d;
960 if (dmax < d) {
961 dmax = d;
962 imax = j;
963 }
964 }
965 if (dmax < m_loopThreshold * dsum) {
966 // If a true loop, close it.
967 contour.push_back(contour[0]);
968 } else {
969 circle = false;
970 if (imax > 0) {
971 // Shift the points to make a line.
972 std::rotate(contour.begin(), contour.begin() + imax, contour.end());
973 }
974 }
975}
976
977}
Abstract base class for components.
virtual bool GetBoundingBox(double &xmin, double &ymin, double &zmin, double &xmax, double &ymax, double &zmax)
Get the bounding box coordinates.
void EnableSignalCalculation(const bool on=true)
Switch on/off calculation of induced currents (default: enabled).
Definition: DriftLineRKF.hh:34
bool DriftNegativeIon(const double x0, const double y0, const double z0, const double t0)
bool DriftIon(const double x0, const double y0, const double z0, const double t0)
Simulate the drift line of an ion with a given starting point.
void SetSensor(Sensor *s)
Set the sensor.
Definition: DriftLineRKF.cc:50
unsigned int GetNumberOfDriftLinePoints() const
Get the number of points of the most recent drift line.
Definition: DriftLineRKF.hh:91
bool DriftElectron(const double x0, const double y0, const double z0, const double t0)
Simulate the drift line of an electron with a given starting point.
bool DriftPositron(const double x0, const double y0, const double z0, const double t0)
void GetDriftLinePoint(const unsigned int i, double &x, double &y, double &z, double &t) const
Get the coordinates and time of a point along the most recent drift line.
void GetEndPoint(double &x, double &y, double &z, double &t, int &st) const
Get the end point and status flag of the most recent drift line.
void SetMaximumStepSize(const double ms)
Set the maximum step size that is allowed. By default, there is no limit.
Definition: DriftLineRKF.cc:67
void AddComponent(ComponentBase *comp)
Add a component.
Definition: Sensor.cc:301
bool GetArea(double &xmin, double &ymin, double &zmin, double &xmax, double &ymax, double &zmax)
Return the current user area.
Definition: Sensor.cc:232
Base class for visualization classes.
Definition: ViewBase.hh:10
std::string m_className
Definition: ViewBase.hh:28
bool m_hasExternalCanvas
Definition: ViewBase.hh:35
TCanvas * m_canvas
Definition: ViewBase.hh:34
void SetArea()
Set the viewing area based on the bounding box of the sensor/component.
void PlotIsochrons(const double tstep, const std::vector< std::array< double, 3 > > &points, const bool rev=false, const bool colour=false, const bool markers=false, const bool plotDriftLines=true)
void SetConnectionThreshold(const double thr)
void SetPlane(const double fx, const double fy, const double fz, const double x0, const double y0, const double z0)
void SetAspectRatioSwitch(const double ar)
ViewIsochrons()
Constructor.
void SetDefaultProjection()
Set the default viewing plane ( - at ).
void SetLoopThreshold(const double thr)
void SetSensor(Sensor *s)
Set the sensor.
void SetComponent(ComponentBase *c)
Set the component.
void Rotate(const double angle)
Rotate the viewing plane (angle in radian).
DoubleAc cos(const DoubleAc &f)
Definition: DoubleAc.cpp:432
DoubleAc fabs(const DoubleAc &f)
Definition: DoubleAc.h:615
DoubleAc sin(const DoubleAc &f)
Definition: DoubleAc.cpp:384
DoubleAc sqrt(const DoubleAc &f)
Definition: DoubleAc.cpp:314