DriftLineRKF.hh

Go to the documentation of this file.

#ifndef G_DRIFTLINE_RKF_H
#define G_DRIFTLINE_RKF_H
 
#include <array>
#include <string>
#include <utility>
#include <vector>
 
#include "Garfield/GarfieldConstants.hh"
 
namespace Garfield {
 
class Sensor;
class ViewDrift;
class Medium;
class DriftLineRKF {
 public:
  DriftLineRKF() : DriftLineRKF(nullptr) {}
  DriftLineRKF(Sensor* sensor);
  ~DriftLineRKF() {}

  void SetSensor(Sensor* s);

  void EnablePlotting(ViewDrift* view);
  void DisablePlotting();

  void EnableSignalCalculation(const bool on = true) { m_doSignal = on; }
  void SetSignalAveragingOrder(const unsigned int navg) { m_navg = navg; }
  void UseWeightingPotential(const bool on = true) {
    m_useWeightingPotential = on;
  }

  void SetIntegrationAccuracy(const double eps);
  void SetMaximumStepSize(const double ms);
  void SetMaximumStepSize();
  void UnsetMaximumStepSize() { m_useStepSizeLimit = false; }
  void RejectKinks(const bool on = true) { m_rejectKinks = on; }

  void SetElectronSignalScalingFactor(const double scale) { m_scaleE = scale; }
  void SetHoleSignalScalingFactor(const double scale) { m_scaleH = scale; }
  void SetIonSignalScalingFactor(const double scale) { m_scaleI = scale; }

  void EnableAvalanche(const bool on = true) { m_doAvalanche = on; }
  void EnableIonTail(const bool on = true) {
    m_doIonTail = on;
    m_doIonTailAuto = false;
  }

  void EnableNegativeIonTail(const bool on = true) { m_doNegativeIonTail = on; }
  void SetGainFluctuationsFixed(const double gain = -1.);
  void SetGainFluctuationsPolya(const double theta, const double mean = -1.,
                                const bool quiet = false);

  void EnableTownsendMap(const bool on = true) { m_useTownsendMap = on; }
  void EnableVelocityMap(const bool on = true) { m_useVelocityMap = on; }

  bool DriftElectron(const double x, const double y, const double z,
                     const double t, const size_t w = 1);
  bool DriftHole(const double x, const double y, const double z,
                 const double t, const size_t w = 1);
  bool DriftIon(const double x, const double y, const double z,
                const double t, const size_t w = 1);
  bool DriftPositron(const double x, const double y, const double z,
                     const double t, const size_t w = 1);
  bool DriftNegativeIon(const double x, const double y, const double z,
                        const double t, const size_t w = 1);

  void PrintDriftLine() const;
  void GetEndPoint(double& x, double& y, double& z, double& t, int& st) const;
  size_t GetNumberOfDriftLinePoints() const { return m_x.size(); }
  void GetDriftLinePoint(const size_t i, double& x, double& y, double& z,
                         double& t) const;

  double GetArrivalTimeSpread(const double eps = 1.e-4) const;
  double GetGain(const double eps = 1.e-4) const;
  double GetLoss(const double eps = 1.e-4) const;
  double GetDriftTime() const {
    return m_t.empty() ? 0. : m_t.back() - m_t.front();
  }

  double GetPathLength() const;

  void GetAvalancheSize(double& ne, double& ni) const {
    ne = m_nE;
    ni = m_nI;
  }

  std::pair<double, double> GetAvalancheSize() const {
    return std::make_pair(m_nE, m_nI);
  }

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

  void EnableDebugging(const bool on = true) { m_debug = on; }
 
 private:
  std::string m_className = "DriftLineRKF";
 
  // Pointer to sensor.
  Sensor* m_sensor = nullptr;
 
  // Particle type of the most recent drift line.
  Particle m_particle = Particle::Electron;
 
  // Maximum allowed step size.
  double m_maxStepSize = 0.;
  // Precision of the stepping algorithm.
  double m_accuracy = 1.e-8;
  // Flag to reject bends > 90 degrees or not.
  bool m_rejectKinks = true;
  // Flag to apply a cut on the maximum allowed step size or not.
  bool m_useStepSizeLimit = false;
 
  // Pointer to the drift viewer.
  ViewDrift* m_view = nullptr;
 
  // Points along the current drift line.
  std::vector<std::array<double, 3> > m_x;
  // Times corresponding to the points along the current drift line.
  std::vector<double> m_t;
  // Status flag of the current drift line.
  int m_status = 0;
 
  // Flag whether to calculate induced signals or not.
  bool m_doSignal = true;
  // Averaging order used when projecting the signal on the time bins.
  unsigned int m_navg = 2;
  // Use weighting potential or weighting field for calculating the signal.
  bool m_useWeightingPotential = true;
 
  // Scaling factor for electron signals.
  double m_scaleE = 1.;
  // Scaling factor for hole signals.
  double m_scaleH = 1.;
  // Scaling factor for ion signals.
  double m_scaleI = 1.;
 
  // Use drift velocity maps?
  bool m_useVelocityMap = false;
  // Use maps for the Townsend coefficient?
  bool m_useTownsendMap = false;
 
  // Flag wether to simulate electron multiplication or not.
  bool m_doAvalanche = true;
  enum class GainFluctuations { None = 0, Polya };
  // Model to be used for randomizing the avalanche size.
  GainFluctuations m_gainFluctuations = GainFluctuations::None;
  // Polya shape parameter.
  double m_theta = 0.;
  // Mean avalanche size (only used if > 1).
  double m_gain = -1.;
 
  // Flag whether to simulate the ion tail or not.
  bool m_doIonTail = true;
  // Simulate the ion tail automatically if the medium has mobility data?
  bool m_doIonTailAuto = true;
  // Flag whether to simulate the negative ion tail or not.
  bool m_doNegativeIonTail = false;
  // Avalanche size.
  double m_nE = 0., m_nI = 0.;
 
  // Debug flag.
  bool m_debug = false;
 
  // Calculate a drift line starting at a given position.
  bool DriftLine(const std::array<double, 3>& x0, const double t0,
                 const Particle particle, std::vector<double>& ts,
                 std::vector<std::array<double, 3> >& xs, int& status) const;
  // Calculate the number of electrons and ions at each point along a
  // drift line.
  bool Avalanche(const Particle particle,
                 const std::vector<std::array<double, 3> >& xs,
                 std::vector<double>& ne, std::vector<double>& ni,
                 std::vector<double>& nn, double& scale) const;
  // Calculate drift lines and induced signals of the positive ions
  // produced in an avalanche.
  bool AddIonTail(const std::vector<double>& te,
                  const std::vector<std::array<double, 3> >& xe,
                  const std::vector<double>& ni, const double scale) const;
  // Calculate drift lines and induced signals of the negative ions
  // produced by electron attachment.
  bool AddNegativeIonTail(const std::vector<double>& te,
                          const std::vector<std::array<double, 3> >& xe,
                          const std::vector<double>& nn,
                          const double scale) const;
  // Compute electric and magnetic field at a given position.
  int GetField(const std::array<double, 3>& x, double& ex, double& ey,
               double& ez, double& bx, double& by, double& bz,
               Medium*& medium) const;
  // Calculate transport parameters for a given point and particle type.
  std::array<double, 3> GetVelocity(const std::array<double, 3>& x,
                                    const Particle particle, int& status) const;
  bool GetDiffusion(const std::array<double, 3>& x, const Particle particle,
                    double& dl, double& dt) const;
  double GetVar(const std::array<double, 3>& x, const Particle particle) const;
  double GetAlpha(const std::array<double, 3>& x,
                  const Particle particle) const;
  double GetEta(const std::array<double, 3>& x, const Particle particle) const;
 
  // Terminate a drift line at the edge of a boundary.
  bool Terminate(const std::array<double, 3>& xx0,
                 const std::array<double, 3>& xx1, const Particle particle,
                 std::vector<double>& ts,
                 std::vector<std::array<double, 3> >& xs) const;
 
  // Drift a particle to a wire
  bool DriftToWire(const double xw, const double yw, const double rw,
                   const Particle particle, std::vector<double>& ts,
                   std::vector<std::array<double, 3> >& xs, int& stat) const;
  // Compute the arrival time spread along a drift line.
  double ComputeSigma(const std::vector<std::array<double, 3> >& x,
                      const Particle particle, const double eps) const;
  // Compute the gain factor along a drift line.
  double ComputeGain(const std::vector<std::array<double, 3> >& x,
                     const Particle particle, const double eps) const;
  // Compute the loss factor along a drift line.
  double ComputeLoss(const std::vector<std::array<double, 3> >& x,
                     const Particle particle, const double eps) const;
  // Integrate the longitudinal diffusion over a step.
  double IntegrateDiffusion(const std::array<double, 3>& xi,
                            const std::array<double, 3>& xe,
                            const Particle particle, const double tol) const;
  // Integrate the Townsend coefficient over a step.
  double IntegrateAlpha(const std::array<double, 3>& xi,
                        const std::array<double, 3>& xe,
                        const Particle particle, const double tol) const;
  // Integrate the attachment coefficient over a step.
  double IntegrateEta(const std::array<double, 3>& xi,
                      const std::array<double, 3>& xe, const Particle particle,
                      const double tol) const;
 
  // Calculate the signal for a given drift line.
  void ComputeSignal(const Particle particle, const double scale,
                     const std::vector<double>& ts,
                     const std::vector<std::array<double, 3> >& xs,
                     const std::vector<double>& ne) const;
 
  // Terminate a field line.
  void Terminate(const std::array<double, 3>& xx0,
                 const std::array<double, 3>& xx1,
                 std::vector<std::array<float, 3> >& xs) const;
 
  static double Charge(const Particle particle) {
    if (particle == Particle::Electron || particle == Particle::NegativeIon) {
      return -1.;
    }
    return 1.;
  }
};
}  // namespace Garfield
 
#endif