#include <DriftLineRKF.hh>

Public Member Functions
	DriftLineRKF ()
	Constructor.

	~DriftLineRKF ()
	Destructor.

void	SetSensor (Sensor *s)
	Set the sensor.

void	EnablePlotting (ViewDrift *view)
	Switch on drift line plotting.

void	DisablePlotting ()
	Switch off drift line plotting.

void	EnableSignalCalculation (const bool on=true)
	Switch on/off calculation of induced currents (default: enabled).

void	SetSignalAveragingOrder (const unsigned int navg)

void	SetIntegrationAccuracy (const double eps)
	Set the accuracy of the Runge Kutta Fehlberg drift line integration.

void	SetMaximumStepSize (const double ms)
	Set the maximum step size that is allowed. By default, there is no limit.

void	UnsetMaximumStepSize ()
	Do not apply an upper limit to the step size that is allowed.

void	RejectKinks (const bool on=true)

void	SetElectronSignalScalingFactor (const double scale)
	Set multiplication factor for the signal induced by electrons.

void	SetHoleSignalScalingFactor (const double scale)
	Set multiplication factor for the signal induced by holes.

void	SetIonSignalScalingFactor (const double scale)
	Set multiplication factor for the signal induced by ions.

void	EnableAvalanche (const bool on=true)
	Enable/disable simulation electron multiplication (default: on).

void	EnableIonTail (const bool on=true)
	Enable/disable simulation of the ion tail (default: off).

void	SetGainFluctuationsFixed (const double gain=-1.)
	Do not randomize the avalanche size.

void	SetGainFluctuationsPolya (const double theta, const double mean=-1.)

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	DriftHole (const double x0, const double y0, const double z0, const double t0)
	Simulate the drift line of a hole with a given starting point.

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.

bool	DriftPositron (const double x0, const double y0, const double z0, const double t0)

bool	DriftNegativeIon (const double x0, const double y0, const double z0, const double t0)

void	PrintDriftLine () const
	Print the trajectory of 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.

unsigned int	GetNumberOfDriftLinePoints () const
	Get the number of points of the most recent drift line.

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.

double	GetArrivalTimeSpread (const double eps=1.e-4)

double	GetGain (const double eps=1.e-4)
	Compute the multiplication factor for the current drift line.

double	GetLoss (const double eps=1.e-4)
	Compute the attachment loss factor for the current drift line.

double	GetDriftTime () const

void	GetAvalancheSize (double &ne, double &ni) const

bool	FieldLine (const double xi, const double yi, const double zi, std::vector< std::array< float, 3 > > &xl, const bool electron=true)

void	EnableDebugging ()

void	DisableDebugging ()

Detailed Description

Calculation of drift lines based on macroscopic transport coefficients using Runge-Kutta-Fehlberg integration.

Definition at line 17 of file DriftLineRKF.hh.

Constructor & Destructor Documentation

◆ DriftLineRKF()

Garfield::DriftLineRKF::DriftLineRKF ( )

Constructor.

Definition at line 50 of file DriftLineRKF.cc.

                           {
  m_t.reserve(1000);
  m_x.reserve(1000);
}

◆ ~DriftLineRKF()

Garfield::DriftLineRKF::~DriftLineRKF ( )

inline

Destructor.

Definition at line 22 of file DriftLineRKF.hh.

22{}

Member Function Documentation

◆ DisableDebugging()

void Garfield::DriftLineRKF::DisableDebugging ( )

inline

Definition at line 116 of file DriftLineRKF.hh.

116{ m_debug = false; }

◆ DisablePlotting()

void Garfield::DriftLineRKF::DisablePlotting ( )

Switch off drift line plotting.

Definition at line 90 of file DriftLineRKF.cc.

90{ m_view = nullptr; }

◆ DriftElectron()

bool Garfield::DriftLineRKF::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.

Definition at line 127 of file DriftLineRKF.cc.

                                                                   {
  m_particle = Particle::Electron;
  if (!DriftLine(x0, y0, z0, t0, Particle::Electron, m_t, m_x, m_status)) {
    return false;
  }
  std::vector<double> ne(m_t.size(), 1.);
  std::vector<double> ni(m_t.size(), 0.);
  double scale = 1.;
  if (m_doAvalanche) Avalanche(Particle::Electron, m_x, ne, ni, scale);
  if (m_doSignal) {
    ComputeSignal(Particle::Electron, scale * m_scaleE, m_t, m_x, ne);
  }
  if (m_doAvalanche && m_doIonTail) AddIonTail(m_t, m_x, ni, scale);
  return true;
}

◆ DriftHole()

bool Garfield::DriftLineRKF::DriftHole	(	const double	x0,
		const double	y0,
		const double	z0,
		const double	t0
	)

Simulate the drift line of a hole with a given starting point.

Definition at line 191 of file DriftLineRKF.cc.

                                              {
  m_particle = Particle::Hole;
  if (!DriftLine(x0, y0, z0, t0, Particle::Hole, m_t, m_x, m_status)) {
    return false;
  }
  if (m_doSignal) ComputeSignal(Particle::Hole, m_scaleH, m_t, m_x, {});
  return true;
}

◆ DriftIon()

bool Garfield::DriftLineRKF::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.

Definition at line 201 of file DriftLineRKF.cc.

                                             {
  m_particle = Particle::Ion;
  if (!DriftLine(x0, y0, z0, t0, Particle::Ion, m_t, m_x, m_status)) {
    return false;
  }
  if (m_doSignal) ComputeSignal(Particle::Ion, m_scaleI, m_t, m_x, {});
  return true;
}

◆ DriftNegativeIon()

bool Garfield::DriftLineRKF::DriftNegativeIon	(	const double	x0,
		const double	y0,
		const double	z0,
		const double	t0
	)

Simulate the drift line of an ion with a given starting point, assuming that it has negative charge.

Definition at line 211 of file DriftLineRKF.cc.

                                                                      {
 
  m_particle = Particle::NegativeIon;
  if (!DriftLine(x0, y0, z0, t0, Particle::NegativeIon, m_t, m_x, m_status)) {
    return false;
  }
  if (m_doSignal) {
    ComputeSignal(Particle::NegativeIon, m_scaleI, m_t, m_x, {});
  }
  return true;
}

◆ DriftPositron()

bool Garfield::DriftLineRKF::DriftPositron	(	const double	x0,
		const double	y0,
		const double	z0,
		const double	t0
	)

Simulate the drift line of an electron with a given starting point, assuming that it has positive charge.

Definition at line 178 of file DriftLineRKF.cc.

                                                                   {
 
  m_particle = Particle::Positron;
  if (!DriftLine(x0, y0, z0, t0, Particle::Positron, m_t, m_x, m_status)) {
    return false;
  }
  if (m_doSignal) {
    ComputeSignal(Particle::Positron, m_scaleE, m_t, m_x, {});
  }
  return true;
}

◆ EnableAvalanche()

void Garfield::DriftLineRKF::EnableAvalanche ( const bool on = true )

inline

Enable/disable simulation electron multiplication (default: on).

Definition at line 58 of file DriftLineRKF.hh.

58{ m_doAvalanche = on; }

◆ EnableDebugging()

void Garfield::DriftLineRKF::EnableDebugging ( )

inline

Definition at line 115 of file DriftLineRKF.hh.

115{ m_debug = true; }

◆ EnableIonTail()

void Garfield::DriftLineRKF::EnableIonTail ( const bool on = true )

inline

Enable/disable simulation of the ion tail (default: off).

Definition at line 60 of file DriftLineRKF.hh.

60{ m_doIonTail = on; }

◆ EnablePlotting()

void Garfield::DriftLineRKF::EnablePlotting ( ViewDrift * view )

Switch on drift line plotting.

Definition at line 82 of file DriftLineRKF.cc.

                                                 {
  if (!view) {
    std::cerr << m_className << "::EnablePlotting: Null pointer.\n";
    return;
  }
  m_view = view;
}

◆ EnableSignalCalculation()

void Garfield::DriftLineRKF::EnableSignalCalculation ( const bool on = true )

inline

Switch on/off calculation of induced currents (default: enabled).

Definition at line 33 of file DriftLineRKF.hh.

33{ m_doSignal = on; }

◆ FieldLine()

bool Garfield::DriftLineRKF::FieldLine	(	const double	xi,
		const double	yi,
		const double	zi,
		std::vector< std::array< float, 3 > > &	xl,
		const bool	electron = `true`
	)

Compute an electric field line.

Parameters

xi,yi,zi	starting point
xl	points along the field line
electron	flag to set the direction in which to follow the field

Definition at line 1603 of file DriftLineRKF.cc.

                                                  {
 
  xl.clear();
 
  // Check if the sensor is defined.
  if (!m_sensor) {
    std::cerr << m_className << "::FieldLine: Sensor is not defined.\n";
    return false;
  }
 
  // Get the sensor's bounding box.
  double xmin = 0., xmax = 0.;
  double ymin = 0., ymax = 0.;
  double zmin = 0., zmax = 0.;
  bool bbox = m_sensor->GetArea(xmin, ymin, zmin, xmax, ymax, zmax);
 
  // Make sure the initial position is at a valid location.
  double ex = 0., ey = 0., ez = 0.;
  Medium* medium = nullptr;
  int stat = 0;
  m_sensor->ElectricField(xi, yi, zi, ex, ey, ez, medium, stat);
  if (!medium || stat != 0) {
    std::cerr << m_className << "::FieldLine:\n"
              << "    No valid field at initial position.\n";
    return false;
  }
  Vec x0 = {xi, yi, zi};
  Vec f0 = {ex, ey, ez};
  if (electron) for (auto& f : f0) f *= -1; 
 
  // Set the numerical constants for the RKF integration.
  constexpr double c10 = 214. / 891.;
  constexpr double c11 = 1. / 33.;
  constexpr double c12 = 650. / 891.;
  constexpr double c20 = 533. / 2106.;
  constexpr double c22 = 800. / 1053.;
  constexpr double c23 = -1. / 78.;
 
  constexpr double b10 = 1. / 4.;
  constexpr double b20 = -189. / 800.;
  constexpr double b21 = 729. / 800.;
  constexpr double b30 = 214. / 891.;
  constexpr double b31 = 1. / 33.;
  constexpr double b32 = 650. / 891.;
 
  const double fmag0 = Mag(f0);
  if (fmag0 < Small) {
    std::cerr << m_className << "::FieldLine:\n"
              << "    Zero field at initial position.\n";
    return false;
  }
 
  // Initialise time step and previous time step.
  double h = m_accuracy / fmag0;
  double hprev = h;
 
  // Set the initial point.
  xl.push_back({float(x0[0]), float(x0[1]), float(x0[2])});
 
  int initCycle = 3;
  bool ok = true;
  while (ok) {
    // Get the field at the first probe point.
    Vec x1 = x0;
    for (unsigned int i = 0; i < 3; ++i) {
      x1[i] += h * b10 * f0[i];
    }
    m_sensor->ElectricField(x1[0], x1[1], x1[2], ex, ey, ez, medium, stat);
    if (stat != 0) {
      if (m_debug) {
        std::cout << m_className << "::FieldLine: Point 1 outside.\n";
      }
      Terminate(x0, x1, xl);
      return true;
    }
    Vec f1 = {ex, ey, ez};
    if (electron) for (auto& f : f1) f *= -1;
    // Get the field at the second probe point.
    Vec x2 = x0;
    for (unsigned int i = 0; i < 3; ++i) {
      x2[i] += h * (b20 * f0[i] + b21 * f1[i]);
    }
    m_sensor->ElectricField(x2[0], x2[1], x2[2], ex, ey, ez, medium, stat);
    if (stat != 0) {
      if (m_debug) {
        std::cout << m_className << "::FieldLine: Point 2 outside.\n";
      }
      Terminate(x0, x2, xl);
      return true;
    }
    Vec f2 = {ex, ey, ez};
    if (electron) for (auto& f : f2) f *= -1;
    // Get the field at the third probe point.
    Vec x3 = x0;
    for (unsigned int i = 0; i < 3; ++i) {
      x3[i] += h * (b30 * f0[i] + b31 * f1[i] + b32 * f2[i]);
    }
    m_sensor->ElectricField(x3[0], x3[1], x3[2], ex, ey, ez, medium, stat);
    if (stat != 0) {
      if (m_debug) {
        std::cout << m_className << "::FieldLine: Point 3 outside.\n";
      }
      Terminate(x0, x3, xl);
      return true;
    }
    Vec f3 = {ex, ey, ez};
    if (electron) for (auto& f : f3) f *= -1;
    // Check if we crossed a wire.
    double xw = 0., yw = 0., zw = 0., rw = 0.;
    if (m_sensor->IsWireCrossed(x0[0], x0[1], x0[2], 
                                x1[0], x1[1], x1[2], xw, yw, zw, false, rw) ||
        m_sensor->IsWireCrossed(x0[0], x0[1], x0[2], 
                                x2[0], x2[1], x2[2], xw, yw, zw, false, rw) ||
        m_sensor->IsWireCrossed(x0[0], x0[1], x0[2], 
                                x3[0], x3[1], x3[2], xw, yw, zw, false, rw)) {
      // TODO!
      xl.push_back({float(xw), float(yw), float(zw)});
      return true;
      // Drift to wire.
      if (h > Small) {
        h *= 0.5;
        continue;
      } else {
        std::cerr << m_className << "::FieldLine: Step size too small. Stop.\n";
        return false;
      }
    }
    // Calculate the correction terms.
    Vec phi1 = {0., 0., 0.};
    Vec phi2 = {0., 0., 0.};
    for (unsigned int i = 0; i < 3; ++i) {
      phi1[i] = c10 * f0[i] + c11 * f1[i] + c12 * f2[i];
      phi2[i] = c20 * f0[i] + c22 * f2[i] + c23 * f3[i];
    }
    // Check if the step length is valid.
    const double phi1mag = Mag(phi1);
    if (phi1mag < Small) {
      std::cerr << m_className << "::FieldLine: Step has zero length. Stop.\n";
      break;
    } else if (m_useStepSizeLimit && h * phi1mag > m_maxStepSize) {
      if (m_debug) {
        std::cout << m_className << "::FieldLine: Step is considered too long. "
                  << "H is reduced.\n";
      }
      h = 0.5 * m_maxStepSize / phi1mag;
      continue;
    } else if (bbox) {
      // Don't allow h to become too large.
      if (h * fabs(phi1[0]) > 0.1 * fabs(xmax - xmin) ||
          h * fabs(phi1[1]) > 0.1 * fabs(ymax - ymin)) {
        h *= 0.5;
        if (m_debug) {
          std::cout << m_className << "::FieldLine: Step is considered too long. "
                    << "H is divided by two.\n";
        }
        continue;
      }
    }
    if (m_debug) std::cout << m_className << "::FieldLine: Step size ok.\n";
    // Update the position.
    for (unsigned int i = 0; i < 3; ++i) x0[i] += h * phi1[i];
    // Check the new position.
    m_sensor->ElectricField(x0[0], x0[1], x0[2], ex, ey, ez, medium, stat);
    if (stat != 0) {
      // The new position is not inside a valid drift medium.
      // Terminate the drift line.
      if (m_debug) {
        std::cout << m_className << "::FieldLine: Point outside. Terminate.\n";
      }
      std::array<double, 3> xp = {xl.back()[0], xl.back()[1], xl.back()[2]};
      Terminate(xp, x0, xl);
      return true;
    }
    // Add the new point to the drift line.
    xl.push_back({float(x0[0]), float(x0[1]), float(x0[2])});
    // Adjust the step size according to the accuracy of the two estimates.
    hprev = h;
    const double dphi = fabs(phi1[0] - phi2[0]) + fabs(phi1[1] - phi2[1]) +
                        fabs(phi1[2] - phi2[2]);
    if (dphi > 0) {
      h = sqrt(h * m_accuracy / dphi);
      if (m_debug) {
        std::cout << m_className << "::FieldLine: Adapting H to " << h << ".\n";
      }
    } else {
      h *= 2;
      if (m_debug) {
        std::cout << m_className << "::FieldLine: H increased by factor two.\n";
      }
    }
    // Make sure that H is different from zero; this should always be ok.
    if (h < Small) {
      std::cerr << m_className << "::FieldLine: Step size is zero. Stop.\n";
      return false;
    }
    // Check the initial step size.
    if (initCycle > 0 && h < 0.2 * hprev) {
      if (m_debug) {
        std::cout << m_className << "::FieldLine: Reinitialise step size.\n";
      }
      --initCycle;
      x0 = {xi, yi, zi};
      xl.clear();
      xl.push_back({float(xi), float(yi), float(zi)});
      continue;
    }
    initCycle = 0;
    // Don't allow H to grow too quickly
    if (h > 10 * hprev) {
      h = 10 * hprev;
      if (m_debug) {
        std::cout << m_className << "::FieldLine: H restricted to 10 times "
                  << "the previous value.\n";
      }
    }
    // Stop in case H tends to become too small.
    if (h * (fabs(phi1[0]) + fabs(phi1[1]) + fabs(phi1[2])) < m_accuracy) {
      std::cerr << m_className << "::FieldLine: Step size has become smaller "
                << "than int. accuracy. Stop.\n";
      return false;
    }
    // Update the field.
    f0 = f3;
  }
  return true;
}

Referenced by Garfield::ViewField::PlotFieldLines().

◆ GetArrivalTimeSpread()

double Garfield::DriftLineRKF::GetArrivalTimeSpread ( const double eps = 1.e-4 )

Compute the sigma of the arrival time distribution for the current drift line by integrating the longitudinal diffusion coefficient.

Definition at line 680 of file DriftLineRKF.cc.

                                                          {
 
  // -----------------------------------------------------------------------
  //    DLCDF1 - Routine returning the integrated diffusion coefficient of
  //             the current drift line. The routine uses an adaptive
  //             Simpson integration.
  // -----------------------------------------------------------------------
 
  const unsigned int nPoints = m_x.size();
  // Return straight away if there is only one point.
  if (nPoints < 2) return 0.;
  const Particle particle = m_particle;
 
  // First get a rough estimate of the result.
  double crude = 0.;
  double varPrev = 0.;
  for (unsigned int i = 0; i < nPoints; ++i) {
    // Get the drift velocity and diffusion coefficients at this step.
    double ex = 0., ey = 0., ez = 0.;
    double bx = 0., by = 0., bz = 0.;
    Medium* medium = nullptr;
    if (GetField(m_x[i], ex, ey, ez, bx, by, bz, medium) != 0) {
      std::cerr << m_className << "::GetArrivalTimeSpread:\n"
                << "    Invalid drift line point " << i << ".\n";
      continue;
    }
    Vec v;
    if (!GetVelocity(ex, ey, ez, bx, by, bz, medium, particle, v)) {
      std::cerr << m_className << "::GetArrivalTimeSpread:\n"
              << "    Cannot retrieve drift velocity at point " << i << ".\n";
      continue;
    }
    const double speed = Mag(v);
    if (speed < Small) {
      std::cerr << m_className << "::GetArrivalTimeSpread:\n"
                << "    Zero drift velocity at point " << i << ".\n";
      continue;
    }
    double dl = 0., dt = 0.;
    if (!GetDiffusion(ex, ey, ez, bx, by, bz, medium, particle, dl, dt)) {
      std::cerr << m_className << "::GetArrivalTimeSpread:\n"
              << "    Cannot retrieve diffusion at point " << i << ".\n";
      continue;
    }
    const double sigma = dl / speed;
    const double var = sigma * sigma;
    if (i > 0) {
      const auto& x = m_x[i];
      const auto& xPrev = m_x[i - 1];
      const double d = Mag(x[0] - xPrev[0], x[1] - xPrev[1], x[2] - xPrev[2]);
      crude += 0.5 * d * (var + varPrev);
    }
    varPrev = var;
  }
  crude = sqrt(crude);
 
  const double tol = eps * crude; 
  double sum = 0.;
  for (unsigned int i = 0; i < nPoints - 1; ++i) {
    sum += IntegrateDiffusion(m_x[i], m_x[i + 1], particle, tol);
  }
  return sqrt(sum);
}

◆ GetAvalancheSize()

void Garfield::DriftLineRKF::GetAvalancheSize	(	double &	ne,
		double &	ni
	)		const

inline

Definition at line 104 of file DriftLineRKF.hh.

104{ ne = m_nE; ni = m_nI; }

◆ GetDriftLinePoint()

void Garfield::DriftLineRKF::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.

Definition at line 1234 of file DriftLineRKF.cc.

                                                                 {
  if (i >= m_x.size()) {
    std::cerr << m_className << "::GetDriftLinePoint: Index out of range.\n";
    return;
  }
 
  const auto& p = m_x[i];
  x = p[0];
  y = p[1];
  z = p[2];
  t = m_t[i];
}

◆ GetDriftTime()

double Garfield::DriftLineRKF::GetDriftTime ( ) const

inline

Definition at line 102 of file DriftLineRKF.hh.

102{ return m_t.empty() ? 0. : m_t.back(); }

◆ GetEndPoint()

void Garfield::DriftLineRKF::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.

Definition at line 1219 of file DriftLineRKF.cc.

                                                {
  if (m_x.empty()) {
    x = y = z = t = 0.;
    stat = m_status;
    return;
  }
  const auto& p = m_x.back();
  x = p[0];
  y = p[1];
  z = p[2];
  t = m_t.back();
  stat = m_status;
}

◆ GetGain()

double Garfield::DriftLineRKF::GetGain ( const double eps = 1.e-4 )

Compute the multiplication factor for the current drift line.

Definition at line 744 of file DriftLineRKF.cc.

                                             {
 
  // -----------------------------------------------------------------------
  //    DLCTW1 - Routine returning the multiplication factor for the current
  //             drift line. The routine uses an adaptive Simpson style
  //             integration.
  // -----------------------------------------------------------------------
 
  const unsigned int nPoints = m_x.size();
  // Return straight away if there is only one point.
  if (nPoints < 2) return 1.;
  const Particle particle = m_particle;
  if (particle == Particle::Ion) return 1.;
  if (m_status == StatusCalculationAbandoned) return 1.;
 
  // First get a rough estimate of the result.
  double crude = 0.;
  double alphaPrev = 0.;
  for (unsigned int i = 0; i < nPoints; ++i) {
    // Get the Townsend coefficient at this step.
    double ex = 0., ey = 0., ez = 0.;
    double bx = 0., by = 0., bz = 0.;
    Medium* medium = nullptr;
    if (GetField(m_x[i], ex, ey, ez, bx, by, bz, medium) != 0) {
      std::cerr << m_className << "::GetGain:\n"
                << "    Invalid drift line point " << i << ".\n";
      continue;
    }
    double alpha = 0.;
    if (!GetAlpha(ex, ey, ez, bx, by, bz, medium, particle, alpha)) {
      std::cerr << m_className << "::GetGain:\n"
                << "    Cannot retrieve alpha at point " << i << ".\n";
      continue;
    }
    if (i > 0) {
      const auto& x = m_x[i];
      const auto& xPrev = m_x[i - 1];
      const double d = Mag(x[0] - xPrev[0], x[1] - xPrev[1], x[2] - xPrev[2]);
      crude += 0.5 * d * (alpha + alphaPrev);
    }
    alphaPrev = alpha;
  }
  // Stop if the rough estimate is negligibly small.
  if (crude < Small) return 1.;
 
  // Calculate the integration tolerance based on the rough estimate.
  const double tol = eps * crude;
  double sum = 0.;
  for (unsigned int i = 0; i < nPoints - 1; ++i) {
    sum += IntegrateAlpha(m_x[i], m_x[i + 1], particle, tol);
  }
  return exp(sum);
}

◆ GetLoss()

double Garfield::DriftLineRKF::GetLoss ( const double eps = 1.e-4 )

Compute the attachment loss factor for the current drift line.

Definition at line 798 of file DriftLineRKF.cc.

                                             {
 
  // -----------------------------------------------------------------------
  //    DLCAT1 - Routine returning the attachment losses for the current
  //             drift line. The routine uses an adaptive Simpson style
  //             integration.
  // -----------------------------------------------------------------------
 
  const unsigned int nPoints = m_x.size();
  // Return straight away if there is only one point.
  if (nPoints < 2) return 1.;
  const Particle particle = m_particle;
  if (particle == Particle::Ion) return 1.;
  if (m_status == StatusCalculationAbandoned) return 1.;
 
  // First get a rough estimate of the result.
  double crude = 0.;
  double etaPrev = 0.;
  for (unsigned int i = 0; i < nPoints; ++i) {
    // Get the Townsend coefficient at this step.
    double ex = 0., ey = 0., ez = 0.;
    double bx = 0., by = 0., bz = 0.;
    Medium* medium = nullptr;
    if (GetField(m_x[i], ex, ey, ez, bx, by, bz, medium) != 0) {
      std::cerr << m_className << "::GetLoss:\n"
                << "    Invalid drift line point " << i << ".\n";
      continue;
    }
    double eta = 0.;
    if (!GetEta(ex, ey, ez, bx, by, bz, medium, particle, eta)) {
      std::cerr << m_className << "::GetLoss:\n"
                << "    Cannot retrieve eta at point " << i << ".\n";
      continue;
    }
    if (i > 0) {
      const auto& x = m_x[i];
      const auto& xPrev = m_x[i - 1];
      const double d = Mag(x[0] - xPrev[0], x[1] - xPrev[1], x[2] - xPrev[2]);
      crude += 0.5 * d * (eta + etaPrev);
    }
    etaPrev = eta;
  }
 
  // Calculate the integration tolerance based on the rough estimate.
  const double tol = eps * crude;
  double sum = 0.;
  for (unsigned int i = 0; i < nPoints - 1; ++i) {
    sum += IntegrateEta(m_x[i], m_x[i + 1], particle, tol);
  }
  return exp(-sum);
}

◆ GetNumberOfDriftLinePoints()

unsigned int Garfield::DriftLineRKF::GetNumberOfDriftLinePoints ( ) const

inline

Get the number of points of the most recent drift line.

Definition at line 90 of file DriftLineRKF.hh.

90{ return m_x.size(); }

◆ PrintDriftLine()

void Garfield::DriftLineRKF::PrintDriftLine ( ) const

Print the trajectory of the most recent drift line.

Definition at line 1192 of file DriftLineRKF.cc.

                                        {
 
  std::cout << m_className << "::PrintDriftLine:\n";
  if (m_x.empty()) {
    std::cout << "    No drift line present.\n";
    return;
  }
  if (m_particle == Particle::Electron) {
    std::cout << "    Particle: electron\n";
  } else if (m_particle == Particle::Ion) {
    std::cout << "    Particle: ion\n";
  } else if (m_particle == Particle::Hole) {
    std::cout << "    Particle: hole\n";
  } else {
    std::cout << "    Particle: unknown\n";
  }
  std::cout << "    Status: " << m_status << "\n"
            << "  Step       time [ns]        "
            << "x [cm]          y [cm]          z [cm]\n";
  const unsigned int nPoints = m_x.size();
  for (unsigned int i = 0; i < nPoints; ++i) {
    std::printf("%6d %15.7f %15.7f %15.7f %15.7f\n", 
                i, m_t[i], m_x[i][0], m_x[i][1], m_x[i][2]);
  }
 
}

◆ RejectKinks()

void Garfield::DriftLineRKF::RejectKinks ( const bool on = true )

inline

Request (or not) the drift line calculation to be aborted if the drift line makes a bend sharper than 90 degrees.

Definition at line 48 of file DriftLineRKF.hh.

48{ m_rejectKinks = on; }

◆ SetElectronSignalScalingFactor()

void Garfield::DriftLineRKF::SetElectronSignalScalingFactor ( const double scale )

inline

Set multiplication factor for the signal induced by electrons.

Definition at line 51 of file DriftLineRKF.hh.

51{ m_scaleE = scale; }

◆ SetGainFluctuationsFixed()

void Garfield::DriftLineRKF::SetGainFluctuationsFixed ( const double gain = -1. )

Do not randomize the avalanche size.

Definition at line 92 of file DriftLineRKF.cc.

                                                             {
 
  if (gain > 1.) {
    std::cout << m_className << "::SetGainFluctuationsFixed: "
              << "Avalanche size set to " << gain << ".\n";
  } else {
    std::cout << m_className << "::SetGainFluctuationsFixed:\n    "
              << "Avalanche size will be given by "
              << "the integrated Townsend coefficient.\n";
  }
  m_gain = gain; 
  m_gainFluctuations = GainFluctuations::None;
}

◆ SetGainFluctuationsPolya()

void Garfield::DriftLineRKF::SetGainFluctuationsPolya	(	const double	theta,
		const double	mean = `-1.`
	)

Sample the avalanche size from a Polya distribution with shape parameter theta.

Definition at line 106 of file DriftLineRKF.cc.

                                                               {
 
  if (theta < 0.) {
    std::cerr << m_className << "::SetGainFluctuationsPolya: "
              << "Shape parameter must be >= 0.\n";
    return;
  }  
  if (mean > 1.) {
    std::cout << m_className << "::SetGainFluctuationsPolya: "
              << "Mean avalanche size set to " << mean << ".\n";
  } else {
    std::cout << m_className << "::SetGainFluctuationsPolya:\n    "
              << "Mean avalanche size will be given by "
              << "the integrated Townsend coefficient.\n";
  }
  m_gain = mean;
  m_theta = theta;
  m_gainFluctuations = GainFluctuations::Polya;
} 

◆ SetHoleSignalScalingFactor()

void Garfield::DriftLineRKF::SetHoleSignalScalingFactor ( const double scale )

inline

Set multiplication factor for the signal induced by holes.

Definition at line 53 of file DriftLineRKF.hh.

53{ m_scaleH = scale; }

◆ SetIntegrationAccuracy()

void Garfield::DriftLineRKF::SetIntegrationAccuracy ( const double eps )

Set the accuracy of the Runge Kutta Fehlberg drift line integration.

Definition at line 63 of file DriftLineRKF.cc.

                                                          {
  if (eps > 0.) {
    m_accuracy = eps;
  } else {
    std::cerr << m_className << "::SetIntegrationAccuracy:\n"
              << "    Accuracy must be greater than zero.\n";
  }
}

◆ SetIonSignalScalingFactor()

void Garfield::DriftLineRKF::SetIonSignalScalingFactor ( const double scale )

inline

Set multiplication factor for the signal induced by ions.

Definition at line 55 of file DriftLineRKF.hh.

55{ m_scaleI = scale; }

◆ SetMaximumStepSize()

void Garfield::DriftLineRKF::SetMaximumStepSize ( const double ms )

Set the maximum step size that is allowed. By default, there is no limit.

Definition at line 72 of file DriftLineRKF.cc.

                                                     {
  if (ms > 0.) {
    m_maxStepSize = ms;
    m_useStepSizeLimit = true;
  } else {
    std::cerr << m_className << "::SetMaximumStepSize:\n"
              << "    Step size must be greater than zero.\n";
  }
}

Referenced by Garfield::ViewField::PlotFieldLines().

◆ SetSensor()

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

Set the sensor.

Definition at line 55 of file DriftLineRKF.cc.

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

Referenced by Garfield::ViewField::PlotFieldLines().

◆ SetSignalAveragingOrder()

void Garfield::DriftLineRKF::SetSignalAveragingOrder ( const unsigned int navg )

inline

Set the number of points to be used when averaging the delayed signal vector over a time bin in the Sensor class. The averaging is done with a $2\times navg + 1$ point Newton-Raphson integration. Default: 2.

Definition at line 38 of file DriftLineRKF.hh.

38{ m_navg = navg; }

◆ UnsetMaximumStepSize()

void Garfield::DriftLineRKF::UnsetMaximumStepSize ( )

inline

Do not apply an upper limit to the step size that is allowed.

Definition at line 45 of file DriftLineRKF.hh.

45{ m_useStepSizeLimit = false; }

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

Public Member Functions

Detailed Description

Constructor & Destructor Documentation

◆ DriftLineRKF()

◆ ~DriftLineRKF()

Member Function Documentation

◆ DisableDebugging()

◆ DisablePlotting()

◆ DriftElectron()

◆ DriftHole()

◆ DriftIon()

◆ DriftNegativeIon()

◆ DriftPositron()

◆ EnableAvalanche()

◆ EnableDebugging()

◆ EnableIonTail()

◆ EnablePlotting()

◆ EnableSignalCalculation()

◆ FieldLine()

◆ GetArrivalTimeSpread()

◆ GetAvalancheSize()

◆ GetDriftLinePoint()

◆ GetDriftTime()

◆ GetEndPoint()

◆ GetGain()

◆ GetLoss()

◆ GetNumberOfDriftLinePoints()

◆ PrintDriftLine()

◆ RejectKinks()

◆ SetElectronSignalScalingFactor()

◆ SetGainFluctuationsFixed()

◆ SetGainFluctuationsPolya()

◆ SetHoleSignalScalingFactor()

◆ SetIntegrationAccuracy()

◆ SetIonSignalScalingFactor()

◆ SetMaximumStepSize()

◆ SetSensor()

◆ SetSignalAveragingOrder()

◆ UnsetMaximumStepSize()