19double Interpolate(
const std::vector<double>& y,
20 const std::vector<double>& x,
const double xx) {
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()];
33bool OnLine(
const double x1,
const double y1,
const double x2,
const double y2,
34 const double u,
const double v) {
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);
41 if ((
fabs(x1 - u) <= epsx &&
fabs(y1 - v) <= epsy) ||
42 (
fabs(x2 - u) <= epsx &&
fabs(y2 - v) <= epsy)) {
45 }
else if (
fabs(x1 - x2) <= epsx &&
fabs(y1 - y2) <= epsy) {
49 double xc = 0., yc = 0.;
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);
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);
82 if (
fabs(u - xc) < epsx &&
fabs(v - yc) < epsy) {
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) {
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)) {
98 std::array<std::array<double, 2>, 2> a;
103 const double det = a[0][0] * a[1][1] - a[1][0] * a[0][1];
107 epsx = std::max(epsx, 1.e-10);
108 epsy = std::max(epsy, 1.e-10);
109 if (
fabs(det) < epsx * epsy) {
114 const double aux = a[1][1];
115 a[1][1] = a[0][0] / det;
117 a[1][0] = -a[1][0] / det;
118 a[0][1] = -a[0][1] / det;
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);
123 if (OnLine(x1, y1, x2, y2, xc, yc) && OnLine(u1, v1, u2, v2, xc, yc)) {
141 std::cerr <<
m_className <<
"::SetSensor: Null pointer.\n";
146 m_component =
nullptr;
151 std::cerr <<
m_className <<
"::SetComponent: Null pointer.\n";
160 const double xmax,
const double ymax) {
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";
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;
178 std::cerr <<
m_className <<
"::SetAspectRatioSwitch: Value must be > 0.\n";
186 if (thr < 0. || thr > 1.) {
187 std::cerr <<
m_className <<
"::SetLoopThreshold:\n"
188 <<
" Value must be between 0 and 1.\n";
191 m_loopThreshold = thr;
196 if (thr < 0. || thr > 1.) {
197 std::cerr <<
m_className <<
"::SetConnectionThreshold:\n"
198 <<
" Value must be between 0 and 1.\n";
201 m_connectionThreshold = thr;
205 const std::vector<std::array<double, 3> >& points,
const bool rev,
206 const bool colour,
const bool markers,
const bool plotDriftLines) {
208 if (!m_sensor && !m_component) {
210 <<
" Neither sensor nor component are defined.\n";
214 std::cerr <<
m_className <<
"::PlotIsochrons: Time step must be > 0.\n";
217 if (points.empty()) {
219 <<
" No starting points provided.\n";
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);
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;
238 ComputeDriftLines(tstep, points, driftLines, startPoints, endPoints,
240 const unsigned int nDriftLines = driftLines.size();
241 if (nDriftLines < 2) {
242 std::cerr <<
m_className <<
"::PlotIsochrons: Too few drift lines.\n";
246 std::size_t nContours = 0;
247 for (
const auto& driftLine : driftLines) {
248 nContours = std::max(nContours, driftLine.size());
250 if (nContours == 0) {
251 std::cerr <<
m_className <<
"::PlotIsochrons: No contour lines.\n";
255 std::set<int> allStats;
256 for (
const auto stat : statusCodes) allStats.insert(stat);
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";
269 std::cout <<
" Do not sort contours.\n";
271 if (m_checkCrossings) {
272 std::cout <<
" Check for crossings.\n";
274 std::cout <<
" Do not check for crossings.\n";
277 std::cout <<
" Mark isochron points.\n";
279 std::cout <<
" Draw isochron lines.\n";
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) {
292 const auto col = gStyle->GetColorPalette(
int((ic + 0.99) * colRange));
293 graph.SetLineColor(col);
294 graph.SetMarkerColor(col);
296 for (
int stat : allStats) {
297 std::vector<std::pair<std::array<double, 4>,
unsigned int> > contour;
299 for (
unsigned int k = 0; k < nDriftLines; ++k) {
300 const auto& dl = driftLines[k];
302 if (statusCodes[k] != stat || ic >= dl.size())
continue;
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));
308 if (contour.empty())
continue;
311 if (m_sortContours && !markers) SortContour(contour, circle);
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]);
326 graph.DrawGraph(xp.size(), xp.data(), yp.data(),
"Psame");
330 const double tolx = (m_xmax - m_xmin) * m_connectionThreshold;
331 const double toly = (m_ymax - m_ymin) * m_connectionThreshold;
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) {
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];
352 if (fabs(x1 - x0) > tolx || fabs(y1 - y0) > toly) gap =
true;
355 const auto i0 = contour[i].second;
356 const auto i1 = contour[i + 1].second;
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)) {
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)) {
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],
386 graph.DrawGraph(xp.size(), xp.data(), yp.data(),
"Lsame");
389 graph.DrawGraph(xp.size(), xp.data(), yp.data(),
"Psame");
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");
406 if (!plotDriftLines)
return;
408 graph.SetLineStyle(1);
409 if (m_particle == Particle::Electron) {
410 graph.SetLineColor(kOrange - 3);
412 graph.SetLineColor(kRed + 1);
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]);
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");
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) {
453 const double lx = 0.1 *
fabs(m_xmax - m_xmin);
454 const double ly = 0.1 *
fabs(m_ymax - m_ymin);
457 for (
const auto& point : points) {
458 if (m_particle == Particle::Electron) {
466 drift.
DriftIon(point[0], point[1], point[2], 0.);
473 if (nu < 3)
continue;
475 double xf = 0., yf = 0., zf = 0., tf = 0.;
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) {
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));
494 std::vector<std::array<double, 3> > tab;
496 for (
unsigned int i = 0; i < nSteps; ++i) {
497 const double t = (i + 1) * tstep;
501 std::array<double, 3> step = {Interpolate(xu, tu, t),
502 Interpolate(yu, tu, t),
503 Interpolate(zu, tu, t)};
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));
513 statusCodes.push_back(status);
515 statusCodes.push_back(0);
542void ViewIsochrons::Labels() {
544 strcpy(m_xLabel,
"\0");
547 const double tol = 1.e-4;
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);
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,
" + ");
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);
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);
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,
" + ");
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);
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);
616 strcat(m_xLabel,
" [cm]");
619 strcpy(m_yLabel,
"\0");
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);
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,
" + ");
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);
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);
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,
" + ");
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);
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);
689 strcat(m_yLabel,
" [cm]");
692 strcpy(m_description,
"\0");
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);
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,
" + ");
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);
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);
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,
" + ");
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);
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);
762 sprintf(buf,
" = %g", m_plane[3]);
763 strcat(m_description, buf);
767 <<
" x label: |" << m_xLabel <<
"|\n"
768 <<
" y label: |" << m_yLabel <<
"|\n"
769 <<
" plane: |" << m_description <<
"|\n";
774 const double x0,
const double y0,
const double z0) {
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;
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);
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);
794 m_proj[1][1] = fz / fyz;
795 m_proj[1][2] = -fy / fyz;
801 <<
" Normal vector has zero norm. No new projection set.\n";
808 m_plane[3] = fx * x0 + fy * y0 + fz * z0;
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];
823 for (
int i = 0; i < 3; ++i) {
824 m_proj[0][i] = auxu[i];
825 m_proj[1][i] = auxv[i];
832void ViewIsochrons::SetupCanvas() {
841bool ViewIsochrons::Range() {
843 if (m_hasUserArea)
return true;
848 if (!m_sensor->
GetArea(bbmin[0], bbmin[1], bbmin[2], bbmax[0], bbmax[1],
851 <<
" Sensor area is not defined.\n"
852 <<
" Please set the plot range explicitly (SetArea).\n";
856 if (!m_component->
GetBoundingBox(bbmin[0], bbmin[1], bbmin[2], bbmax[0],
857 bbmax[1], bbmax[2])) {
859 <<
" Bounding box of the component is not defined.\n"
860 <<
" Please set the plot range explicitly (SetArea).\n";
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;
889void ViewIsochrons::SortContour(
890 std::vector<std::pair<std::array<double, 4>,
unsigned int> >& contour,
893 if (contour.size() < 2)
return;
897 for (
const auto& point : contour) {
898 xcog += point.first[0];
899 ycog += point.first[1];
901 const double scale = 1. / contour.size();
908 for (
const auto& point : contour) {
909 const double dx = point.first[0] - xcog;
910 const double dy = point.first[1] - ycog;
915 const double theta = 0.5 * atan2(2 * sxy, sxx - syy);
916 const double ct =
cos(theta);
917 const double st =
sin(theta);
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);
928 if (
fabs(sxx) > m_aspectRatio *
fabs(syy) ||
929 fabs(syy) > m_aspectRatio *
fabs(sxx)) {
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;
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]; }
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;
965 if (dmax < m_loopThreshold * dsum) {
967 contour.push_back(contour[0]);
972 std::rotate(contour.begin(), contour.begin() + imax, contour.end());
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).
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.
unsigned int GetNumberOfDriftLinePoints() const
Get the number of points of the most recent drift line.
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.
void AddComponent(ComponentBase *comp)
Add a component.
bool GetArea(double &xmin, double &ymin, double &zmin, double &xmax, double &ymax, double &zmax)
Return the current user area.
Base class for visualization classes.
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)
DoubleAc fabs(const DoubleAc &f)
DoubleAc sin(const DoubleAc &f)
DoubleAc sqrt(const DoubleAc &f)