Draw equal time contour lines. More...

#include <ViewIsochrons.hh>

Inheritance diagram for Garfield::ViewIsochrons:

Public Member Functions
	ViewIsochrons ()
	Constructor.

	~ViewIsochrons ()=default
	Destructor.

void	SetSensor (Sensor *s)
	Set the sensor.

void	SetComponent (ComponentBase *c)
	Set the component.

void	SetArea (const double xmin, const double ymin, const double xmax, const double ymax)
	Set the viewing area (in local coordinates of the current viewing plane).

void	SetArea ()
	Set the viewing area based on the bounding box of the sensor/component.

void	SetPlane (const double fx, const double fy, const double fz, const double x0, const double y0, const double z0)

void	SetDefaultProjection ()
	Set the default viewing plane ( - at ).

void	Rotate (const double angle)
	Rotate the viewing plane (angle in radian).

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	DriftElectrons (const bool positive=false)

void	DriftIons (const bool negative=false)
	Request ion drift lines with positive (default) or negative charge.

void	EnableSorting (const bool on=true)
	Sort (or not) the points on a contour line (default: sorting is done).

void	CheckCrossings (const bool on=true)

void	SetAspectRatioSwitch (const double ar)

void	SetLoopThreshold (const double thr)

void	SetConnectionThreshold (const double thr)

Public Member Functions inherited from Garfield::ViewBase
	ViewBase ()=delete
	Default constructor.

	ViewBase (const std::string &name)
	Constructor.

virtual	~ViewBase ()
	Destructor.

void	SetCanvas (TCanvas *c)
	Set the canvas to be painted on.

TCanvas *	GetCanvas ()
	Retrieve the canvas.

void	EnableDebugging (const bool on=true)
	Switch on/off debugging output.

Additional Inherited Members
Protected Member Functions inherited from Garfield::ViewBase
std::string	FindUnusedFunctionName (const std::string &s) const

std::string	FindUnusedHistogramName (const std::string &s) const

Protected Attributes inherited from Garfield::ViewBase
std::string	m_className = "ViewBase"

bool	m_debug = false

TCanvas *	m_canvas = nullptr

bool	m_hasExternalCanvas = false

double	m_proj [3][3]

Detailed Description

Draw equal time contour lines.

Definition at line 13 of file ViewIsochrons.hh.

Constructor & Destructor Documentation

◆ ViewIsochrons()

Garfield::ViewIsochrons::ViewIsochrons ( )

Constructor.

Definition at line 135 of file ViewIsochrons.cc.

                             : ViewBase("ViewIsochrons") {
  SetDefaultProjection();
}

◆ ~ViewIsochrons()

Garfield::ViewIsochrons::~ViewIsochrons ( )

default

Destructor.

Member Function Documentation

◆ CheckCrossings()

void Garfield::ViewIsochrons::CheckCrossings ( const bool on = true )

inline

Check (or not) that drift-lines do not cross isochrons (default: check is done).

Definition at line 72 of file ViewIsochrons.hh.

72{ m_checkCrossings = on; }

◆ DriftElectrons()

void Garfield::ViewIsochrons::DriftElectrons ( const bool positive = false )

inline

Request electron drift lines with negative (default) or positive charge.

Definition at line 59 of file ViewIsochrons.hh.

                                                   { 
    m_particle = Particle::Electron;
    m_positive = positive;
  }

◆ DriftIons()

void Garfield::ViewIsochrons::DriftIons ( const bool negative = false )

inline

Request ion drift lines with positive (default) or negative charge.

Definition at line 64 of file ViewIsochrons.hh.

                                              { 
    m_particle = Particle::Ion;
    m_positive = !negative;
  }

◆ EnableSorting()

void Garfield::ViewIsochrons::EnableSorting ( const bool on = true )

inline

Sort (or not) the points on a contour line (default: sorting is done).

Definition at line 69 of file ViewIsochrons.hh.

69{ m_sortContours = on; }

◆ PlotIsochrons()

void Garfield::ViewIsochrons::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`
	)

Draw equal time contour lines.

Parameters

tstep	Time interval.
points	List of starting points.
rev	If true, the drift time is measured from the end points of the drift lines.
colour	Draw contour lines using colours.
markers	Draw markers (as opposed to lines).
plotDriftLines	Draw drift lines together with the isochrons.

Definition at line 204 of file ViewIsochrons.cc.

                                                                      {
 
  if (!m_sensor && !m_component) {
    std::cerr << m_className << "::PlotIsochrons:\n"
              << "    Neither sensor nor component are defined.\n";
    return;
  }
  if (tstep <= 0.) {
    std::cerr << m_className << "::PlotIsochrons: Time step must be > 0.\n";
    return;
  }
  if (points.empty()) {
    std::cerr << m_className << "::PlotIsochrons:\n"
              << "    No starting points provided.\n";
    return;
  }
  SetupCanvas();
  if (!Range()) return;
  auto frame = m_canvas->DrawFrame(m_xmin, m_ymin, m_xmax, m_ymax);
  frame->GetXaxis()->SetTitle(m_xLabel);
  frame->GetYaxis()->SetTitle(m_yLabel);
  m_canvas->Update();
 
  //-----------------------------------------------------------------------
  //   DRFEQT - The main routine (DRFEQT) accumulates equal drift time data
  //   DRFEQP   which is plotted as a set of contours in the entry DRFEQP.
  //-----------------------------------------------------------------------
  std::vector<std::vector<std::array<double, 3> > > driftLines;
  std::vector<std::array<double, 3> > startPoints;  
  std::vector<std::array<double, 3> > endPoints;  
  std::vector<int> statusCodes;
  // Accumulate drift lines.
  ComputeDriftLines(tstep, points, driftLines, startPoints, endPoints, 
                    statusCodes, rev);
  const unsigned int nDriftLines = driftLines.size();
  if (nDriftLines < 2) {
    std::cerr << m_className << "::PlotIsochrons: Too few drift lines.\n";
    return;
  }
  // Keep track of the largest number of contours.
  std::size_t nContours = 0;
  for (const auto& driftLine : driftLines) {
    nContours = std::max(nContours, driftLine.size());
  }
  if (nContours == 0) {
    std::cerr << m_className << "::PlotIsochrons: No contour lines.\n";
    return;
  }
 
  std::set<int> allStats;
  for (const auto stat : statusCodes) allStats.insert(stat);
 
  // DRFEQP
  if (m_debug) {
    std::cout << m_className << "::PlotIsochrons:\n"
              << "    Drawing " << nContours << " contours, "
              << nDriftLines << " drift lines.\n";
    std::printf("    Connection threshold:   %10.3f\n", m_connectionThreshold);
    std::printf("    Aspect ratio threshold: %10.3f\n", m_aspectRatio);
    std::printf("    Loop closing threshold: %10.3f\n", m_loopThreshold);
    if (m_sortContours) {
      std::cout << "    Sort contours.\n";
    } else {
      std::cout << "    Do not sort contours.\n";
    }
    if (m_checkCrossings) {
      std::cout << "    Check for crossings.\n";
    } else {
      std::cout << "    Do not check for crossings.\n";
    }
    if (markers) {
      std::cout << "    Mark isochron points.\n";
    } else {
      std::cout << "    Draw isochron lines.\n";
    }
  }
 
  // Loop over the equal time contours.
  TGraph graph;
  graph.SetLineColor(kGray + 2);
  graph.SetMarkerColor(kGray + 2);
  graph.SetLineStyle(m_lineStyle);
  graph.SetMarkerStyle(m_markerStyle);
  const double colRange = double(gStyle->GetNumberOfColors()) / nContours;
  for (unsigned int ic = 0; ic < nContours; ++ic) {
    if (colour) {
      const auto col = gStyle->GetColorPalette(int((ic + 0.99) * colRange));
      graph.SetLineColor(col);
      graph.SetMarkerColor(col);
    }
    for (int stat : allStats) {
      std::vector<std::pair<std::array<double, 4>, unsigned int> > contour;
      // Loop over the drift lines, picking up the points when OK.
      for (unsigned int k = 0; k < nDriftLines; ++k) {
        const auto& dl = driftLines[k]; 
        // Reject any undesirable combinations.
        if (statusCodes[k] != stat || ic >= dl.size()) continue;
        // Add the point to the contour line.
        std::array<double, 4> point = {dl[ic][0], dl[ic][1], dl[ic][2], 0.};
        contour.push_back(std::make_pair(point, k)); 
      }
      // Skip the plot of this contour if there are no points.
      if (contour.empty()) continue;
      bool circle = false;
      // If requested, sort the points on the contour line.
      if (m_sortContours && !markers) SortContour(contour, circle);
      // Plot this contour.
      if (markers) {
        // Simply mark the contours if this was requested.
        std::vector<double> xp;
        std::vector<double> yp;
        std::vector<double> zp;
        for (const auto& point : contour) {
          const double x = point.first[0];
          const double y = point.first[1];
          const double z = point.first[2];
          xp.push_back(m_proj[0][0] * x + m_proj[1][0] * y + z * m_plane[0]);
          yp.push_back(m_proj[0][1] * x + m_proj[1][1] * y + z * m_plane[1]);
          zp.push_back(m_proj[0][2] * x + m_proj[1][2] * y + z * m_plane[2]);
        }
        graph.DrawGraph(xp.size(), xp.data(), yp.data(), "Psame");
        continue;
      }
      // Regular plotting.
      const double tolx = (m_xmax - m_xmin) * m_connectionThreshold;
      const double toly = (m_ymax - m_ymin) * m_connectionThreshold;
      // Flag to keep track if the segment is interrupted by a drift line
      // or if it is too long.
      bool gap = false;
      // Coordinates to be plotted.
      std::vector<double> xp;
      std::vector<double> yp;
      std::vector<double> zp;
      const unsigned int nP = contour.size();
      for (unsigned int i = 0; i < nP; ++i) { 
        gap = false;
        const auto x0 = contour[i].first[0];
        const auto y0 = contour[i].first[1];
        const auto z0 = contour[i].first[2];
        xp.push_back(m_proj[0][0] * x0 + m_proj[1][0] * y0 + z0 * m_plane[0]);
        yp.push_back(m_proj[0][1] * x0 + m_proj[1][1] * y0 + z0 * m_plane[1]);
        zp.push_back(m_proj[0][2] * x0 + m_proj[1][2] * y0 + z0 * m_plane[2]);
        if (i == nP - 1) break; 
        const auto x1 = contour[i + 1].first[0];
        const auto y1 = contour[i + 1].first[1];
        // Reject contour segments which are long compared with AREA.
        if (fabs(x1 - x0) > tolx || fabs(y1 - y0) > toly) gap = true;
        // Get the indices of the drift lines corresponding 
        // to these two points on the contour line.
        const auto i0 = contour[i].second;
        const auto i1 = contour[i + 1].second; 
        // Set the BREAK flag if it crosses some stored drift line segment.
        if (m_checkCrossings && !gap) {
          for (unsigned int k = 0; k < nDriftLines; ++k) {
            const auto& dl = driftLines[k];
            for (unsigned int jc = 0; jc < dl.size(); ++jc) {
              if ((i0 == k || i1 == k) && (jc == ic || jc + 1 == ic)) {
                continue;
              }
              const auto& p0 = dl[jc];
              const auto& p1 = jc == dl.size() - 1 ? endPoints[k] : dl[jc + 1];
              if (Crossing(p0[0], p0[1], p1[0], p1[1], x0, y0, x1, y1)) {
                gap = true;
                break;
              }
            }
            if (gap) break;
            if ((i0 == k || i1 == k) && ic == 0) continue;
            const auto& p0 = startPoints[k];
            if (Crossing(p0[0], p0[1], dl[0][0], dl[0][1], 
                         x0, y0, x1, y1)) {
              gap = true;
              break;
            }
          }
        }
        // If there has been a break, plot what we have already.
        if (gap) {
          if (xp.size() > 1) {
            // Plot line.
            graph.DrawGraph(xp.size(), xp.data(), yp.data(), "Lsame");
          } else {
            // Plot markers.
            graph.DrawGraph(xp.size(), xp.data(), yp.data(), "Psame");
          }
          xp.clear();
          yp.clear();
          zp.clear();
        }
      }
      // Plot the remainder.
      if (xp.size() > 1) {
        graph.DrawGraph(xp.size(), xp.data(), yp.data(), "Lsame");
      } else if (!xp.empty()) {
        graph.DrawGraph(xp.size(), xp.data(), yp.data(), "Psame");
      }
    }
  }
 
  gPad->Update();
  if (!plotDriftLines) return;
 
  graph.SetLineStyle(1);
  if (m_particle == Particle::Electron) {
    graph.SetLineColor(kOrange - 3);
  } else {
    graph.SetLineColor(kRed + 1);
  }
  for (unsigned int i = 0; i < nDriftLines; ++i) {
    std::vector<double> xp;
    std::vector<double> yp;
    const double x0 = startPoints[i][0];
    const double y0 = startPoints[i][1];
    const double z0 = startPoints[i][2];
    xp.push_back(m_proj[0][0] * x0 + m_proj[1][0] * y0 + z0 * m_plane[0]);
    yp.push_back(m_proj[0][1] * x0 + m_proj[1][1] * y0 + z0 * m_plane[1]);
    for (const auto& point : driftLines[i]) {
      const double x = point[0];
      const double y = point[1];
      const double z = point[2];
      xp.push_back(m_proj[0][0] * x + m_proj[1][0] * y + z * m_plane[0]);
      yp.push_back(m_proj[0][1] * x + m_proj[1][1] * y + z * m_plane[1]);
    }
    const double x1 = endPoints[i][0];
    const double y1 = endPoints[i][1];
    const double z1 = endPoints[i][2];
    xp.push_back(m_proj[0][0] * x1 + m_proj[1][0] * y1 + z1 * m_plane[0]);
    yp.push_back(m_proj[0][1] * x1 + m_proj[1][1] * y1 + z1 * m_plane[1]);
    graph.DrawGraph(xp.size(), xp.data(), yp.data(), "Lsame");
  }
}

◆ Rotate()

void Garfield::ViewIsochrons::Rotate ( const double angle )

Rotate the viewing plane (angle in radian).

Definition at line 814 of file ViewIsochrons.cc.

                                             {
  // Rotate the axes
  double auxu[3], auxv[3];
  const double ctheta = cos(theta);
  const double stheta = sin(theta);
  for (int i = 0; i < 3; ++i) {
    auxu[i] = ctheta * m_proj[0][i] - stheta * m_proj[1][i];
    auxv[i] = stheta * m_proj[0][i] + ctheta * m_proj[1][i];
  }
  for (int i = 0; i < 3; ++i) {
    m_proj[0][i] = auxu[i];
    m_proj[1][i] = auxv[i];
  }
 
  // Make labels to be placed along the axes
  Labels();
}

◆ SetArea() [1/2]

void Garfield::ViewIsochrons::SetArea ( )

inline

Set the viewing area based on the bounding box of the sensor/component.

Definition at line 29 of file ViewIsochrons.hh.

29{ m_hasUserArea = false; }

◆ SetArea() [2/2]

void Garfield::ViewIsochrons::SetArea	(	const double	xmin,
		const double	ymin,
		const double	xmax,
		const double	ymax
	)

Set the viewing area (in local coordinates of the current viewing plane).

Definition at line 159 of file ViewIsochrons.cc.

                                                                  {
  // Check range, assign if non-null.
  if (xmin == xmax || ymin == ymax) {
    std::cerr << m_className << "::SetArea: Null area range is not permitted.\n"
              << "      " << xmin << " < x < " << xmax << "\n"
              << "      " << ymin << " < y < " << ymax << "\n";
    return;
  }
  m_xmin = std::min(xmin, xmax);
  m_ymin = std::min(ymin, ymax);
  m_xmax = std::max(xmin, xmax);
  m_ymax = std::max(ymin, ymax);
  m_hasUserArea = true;
}

◆ SetAspectRatioSwitch()

void Garfield::ViewIsochrons::SetAspectRatioSwitch ( const double ar )

Set the aspect ratio above which an isochron is considered linear (as opposed to circular).

Definition at line 175 of file ViewIsochrons.cc.

                                                        {
 
  if (ar < 0.) {
    std::cerr << m_className << "::SetAspectRatioSwitch: Value must be > 0.\n";
    return;
  }
  m_aspectRatio = ar;
}

◆ SetComponent()

void Garfield::ViewIsochrons::SetComponent ( ComponentBase * c )

Set the component.

Definition at line 149 of file ViewIsochrons.cc.

                                                 {
  if (!c) {
    std::cerr << m_className << "::SetComponent: Null pointer.\n";
    return;
  }
 
  m_component = c;
  m_sensor = nullptr;
}

◆ SetConnectionThreshold()

void Garfield::ViewIsochrons::SetConnectionThreshold ( const double thr )

Fractional distance over which isochron segments are connected (default: 0.2).

Definition at line 194 of file ViewIsochrons.cc.

                                                           {
 
  if (thr < 0. || thr > 1.) {
    std::cerr << m_className << "::SetConnectionThreshold:\n"
              << "    Value must be between 0 and 1.\n";
    return;
  }
  m_connectionThreshold = thr;
}

◆ SetDefaultProjection()

void Garfield::ViewIsochrons::SetDefaultProjection ( )

Set the default viewing plane ( $x$ - $y$ at $z = 0$ ).

Definition at line 520 of file ViewIsochrons.cc.

                                         {
  // Default projection: x-y at z=0
  m_proj[0][0] = 1;
  m_proj[1][0] = 0;
  m_proj[2][0] = 0;
  m_proj[0][1] = 0;
  m_proj[1][1] = 1;
  m_proj[2][1] = 0;
  m_proj[0][2] = 0;
  m_proj[1][2] = 0;
  m_proj[2][2] = 0;
 
  // Plane description
  m_plane[0] = 0;
  m_plane[1] = 0;
  m_plane[2] = 1;
  m_plane[3] = 0;
 
  // Prepare axis labels.
  Labels();
}

Referenced by ViewIsochrons().

◆ SetLoopThreshold()

void Garfield::ViewIsochrons::SetLoopThreshold ( const double thr )

Fractional distance between two points for closing a circular isochron (default: 0.2).

Definition at line 184 of file ViewIsochrons.cc.

                                                     {
 
  if (thr < 0. || thr > 1.) {
    std::cerr << m_className << "::SetLoopThreshold:\n"
              << "    Value must be between 0 and 1.\n";
    return;
  }
  m_loopThreshold = thr;
}

◆ SetPlane()

void Garfield::ViewIsochrons::SetPlane	(	const double	fx,
		const double	fy,
		const double	fz,
		const double	x0,
		const double	y0,
		const double	z0
	)

Set the projection (viewing plane).

Parameters

fx,fy,fz	normal vector
x0,y0,z0	in-plane point

Definition at line 773 of file ViewIsochrons.cc.

                                                                            {
  // Calculate two in-plane vectors for the normal vector
  const double fnorm = sqrt(fx * fx + fy * fy + fz * fz);
  if (fnorm > 0 && fx * fx + fz * fz > 0) {
    const double fxz = sqrt(fx * fx + fz * fz);
    m_proj[0][0] = fz / fxz;
    m_proj[0][1] = 0;
    m_proj[0][2] = -fx / fxz;
    m_proj[1][0] = -fx * fy / (fxz * fnorm);
    m_proj[1][1] = (fx * fx + fz * fz) / (fxz * fnorm);
    m_proj[1][2] = -fy * fz / (fxz * fnorm);
    m_proj[2][0] = x0;
    m_proj[2][1] = y0;
    m_proj[2][2] = z0;
  } else if (fnorm > 0 && fy * fy + fz * fz > 0) {
    const double fyz = sqrt(fy * fy + fz * fz);
    m_proj[0][0] = (fy * fy + fz * fz) / (fyz * fnorm);
    m_proj[0][1] = -fx * fz / (fyz * fnorm);
    m_proj[0][2] = -fy * fz / (fyz * fnorm);
    m_proj[1][0] = 0;
    m_proj[1][1] = fz / fyz;
    m_proj[1][2] = -fy / fyz;
    m_proj[2][0] = x0;
    m_proj[2][1] = y0;
    m_proj[2][2] = z0;
  } else {
    std::cout << m_className << "::SetPlane:\n"
              << "    Normal vector has zero norm. No new projection set.\n";
  }
 
  // Store the plane description
  m_plane[0] = fx;
  m_plane[1] = fy;
  m_plane[2] = fz;
  m_plane[3] = fx * x0 + fy * y0 + fz * z0;
 
  // Make labels to be placed along the axes
  Labels();
}

◆ SetSensor()

void Garfield::ViewIsochrons::SetSensor ( Sensor * s )

Set the sensor.

Definition at line 139 of file ViewIsochrons.cc.

                                       {
  if (!s) {
    std::cerr << m_className << "::SetSensor: Null pointer.\n";
    return;
  }
 
  m_sensor = s;
  m_component = nullptr;
}

The documentation for this class was generated from the following files:

Public Member Functions

Additional Inherited Members

Detailed Description

Constructor & Destructor Documentation

◆ ViewIsochrons()

◆ ~ViewIsochrons()

Member Function Documentation

◆ CheckCrossings()

◆ DriftElectrons()

◆ DriftIons()

◆ EnableSorting()

◆ PlotIsochrons()

◆ Rotate()

◆ SetArea() [1/2]

◆ SetArea() [2/2]

◆ SetAspectRatioSwitch()

◆ SetComponent()

◆ SetConnectionThreshold()

◆ SetDefaultProjection()

◆ SetLoopThreshold()

◆ SetPlane()

◆ SetSensor()