AvalancheMC.hh

Go to the documentation of this file.

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

  void SetSensor(Sensor* s);

  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 DriftNegativeIon(const double x, const double y, const double z,
                        const double t, const size_t w = 1);
  bool AvalancheElectron(const double x, const double y, const double z,
                         const double t, const bool hole = false,
                         const size_t w = 1);
  bool AvalancheHole(const double x, const double y, const double z,
                     const double t, const bool electron = false,
                     const size_t w = 1);
  bool AvalancheElectronHole(const double x, const double y, const double z,
                             const double t, const size_t w = 1);

  void AddElectron(const double x, const double y, const double z,
                   const double t, const size_t w = 1);
  void AddHole(const double x, const double y, const double z, 
               const double t, const size_t w = 1);
  void AddIon(const double x, const double y, const double z, 
              const double t, const size_t w = 1);
  void AddNegativeIon(const double x, const double y, const double z,
                      const double t, const size_t w = 1);
  bool ResumeAvalanche(const bool electron = true, const bool hole = true);
 
  struct Point {
    double x, y, z;
    double t;
  };
 
  struct EndPoint {
    int status;               
    std::vector<Point> path;  
    size_t weight;            
  };
 
  struct Seed {
    Point pt;      
    Particle type; 
    size_t w = 1;  
  };
 
  const std::vector<EndPoint>& GetElectrons() const { return m_electrons; }
  const std::vector<EndPoint>& GetHoles() const { return m_holes; }
  const std::vector<EndPoint>& GetIons() const { return m_ions; }
  const std::vector<EndPoint>& GetNegativeIons() const {
    return m_negativeIons;
  }

  size_t GetNumberOfElectronEndpoints() const { return m_electrons.size(); }
  size_t GetNumberOfIonEndpoints() const { return m_ions.size(); }

  void GetElectronEndpoint(const size_t i, double& x0, double& y0, double& z0,
                           double& t0, double& x1, double& y1, double& z1,
                           double& t1, int& status) const;
  void GetIonEndpoint(const size_t i, double& x0, double& y0, double& z0,
                      double& t0, double& x1, double& y1, double& z1,
                      double& t1, int& status) const;
  void GetNegativeIonEndpoint(const size_t i, double& x0, double& y0,
                              double& z0, double& t0, double& x1, double& y1,
                              double& z1, double& t1, int& status) const;

  void EnablePlotting(ViewDrift* view);
  void DisablePlotting() { m_viewer = nullptr; }

  void EnableDriftLines(const bool on = true) { m_storeDriftLines = on; }

  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 EnableInducedChargeCalculation(const bool on = true) {
    m_doInducedCharge = on;
  }

  void EnableRKFSteps(const bool on = true) { m_doRKF = on; }

  void EnableProjectedPathIntegration(const bool on = true) {
    m_doEquilibration = on;
  }

  void EnableDiffusion() { m_useDiffusion = true; }
  void DisableDiffusion() { m_useDiffusion = false; }

  void EnableAttachment() { m_useAttachment = true; }
  void DisableAttachment() { m_useAttachment = false; }
  void EnableMultiplication(const bool on) { m_useMultiplication = on; }

  void EnableTownsendMap(const bool on = true) { m_useTownsendMap = on; }
  void EnableAttachmentMap(const bool on = true) { m_useAttachmentMap = on; }
  void EnableMobilityMap(const bool on = true) { m_useMobilityMap = on; }
  void EnableVelocityMap(const bool on = true) { m_useVelocityMap = on; }

  void EnableAvalancheSizeLimit(const unsigned int size) { m_sizeCut = size; }
  void DisableAvalancheSizeLimit() { m_sizeCut = 0; }
  int GetAvalancheSizeLimit() const { return m_sizeCut; }

  void SetTimeSteps(const double d = 0.02);
  void SetDistanceSteps(const double d = 0.001);
  void SetCollisionSteps(const unsigned int n = 100);
  void SetStepDistanceFunction(double (*f)(double x, double y, double z));

  void SetTimeWindow(const double t0, const double t1);
  void UnsetTimeWindow() { m_hasTimeWindow = false; }

  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 GetAvalancheSize(unsigned int& ne, unsigned int& ni) const {
    ne = m_nElectrons;
    ni = std::max(m_nIons, m_nHoles);
  }

  std::pair<unsigned int, unsigned int> GetAvalancheSize() const {
    return std::make_pair(m_nElectrons, std::max(m_nIons, m_nHoles));
  }

  void EnableDebugging(const bool on = true) { m_debug = on; }
 
 private:
  std::string m_className = "AvalancheMC";
 
  Sensor* m_sensor = nullptr;
 
  enum class StepModel {
    FixedTime,
    FixedDistance,
    CollisionTime,
    UserDistance
  };

  StepModel m_stepModel = StepModel::CollisionTime;

  double m_tMc = 0.02;
  double m_dMc = 0.001;
  int m_nMc = 100;
  double (*m_fStep)(double x, double y, double z) = nullptr;

  bool m_hasTimeWindow = false;
  double m_tMin = 0.;
  double m_tMax = 0.;

  unsigned int m_sizeCut = 0;

  unsigned int m_nElectrons = 0;
  unsigned int m_nHoles = 0;
  unsigned int m_nIons = 0;
  unsigned int m_nNegativeIons = 0;

  std::vector<EndPoint> m_electrons;
  std::vector<EndPoint> m_holes;
  std::vector<EndPoint> m_ions;
  std::vector<EndPoint> m_negativeIons;
 
  ViewDrift* m_viewer = nullptr;
 
  bool m_storeDriftLines = false;
  bool m_doSignal = true;
  unsigned int m_navg = 1;
  bool m_useWeightingPotential = true;
  bool m_doInducedCharge = false;
  bool m_doEquilibration = true;
  bool m_doRKF = false;
  bool m_useDiffusion = true;
  bool m_useAttachment = true;
  bool m_useMultiplication = true;
  double m_scaleE = 1.;
  double m_scaleH = 1.;
  double m_scaleI = 1.;

  bool m_useTownsendMap = false;
  bool m_useAttachmentMap = false;
  bool m_useMobilityMap = false;
  bool m_useVelocityMap = false;
 
  bool m_debug = false;

  int DriftLine(const Seed& seed, std::vector<Point>& path,
                std::vector<Seed>& secondaries,
                const bool aval, const bool signal) const;
  bool TransportParticles(std::vector<Seed>& stack,
                          const bool withElectrons, const bool withHoles,
                          const bool aval);

  int GetField(const std::array<double, 3>& x, std::array<double, 3>& e,
               std::array<double, 3>& b, Medium*& medium) const;
  double GetMobility(const Particle particle, Medium* medium,
                     const std::array<double, 3>& x) const;
  bool GetVelocity(const Particle particle, Medium* medium,
                   const std::array<double, 3>& x,
                   const std::array<double, 3>& e,
                   const std::array<double, 3>& b,
                   std::array<double, 3>& v) const;
  bool GetDiffusion(const Particle particle, Medium* medium,
                    const std::array<double, 3>& e,
                    const std::array<double, 3>& b, double& dl,
                    double& dt) const;
  double GetAttachment(const Particle particle, Medium* medium,
                       const std::array<double, 3>& x,
                       const std::array<double, 3>& e,
                       const std::array<double, 3>& b) const;
  double GetTownsend(const Particle particle, Medium* medium,
                     const std::array<double, 3>& x,
                     const std::array<double, 3>& e,
                     const std::array<double, 3>& b) const;
  void StepRKF(const Particle particle, const std::array<double, 3>& x0,
               const std::array<double, 3>& v0, const double dt,
               std::array<double, 3>& xf, std::array<double, 3>& vf,
               int& status) const;
  void AddDiffusion(const double step, const double dl, const double dt,
                    std::array<double, 3>& x,
                    const std::array<double, 3>& v) const;
  void Terminate(const std::array<double, 3>& x0, const double t0,
                 std::array<double, 3>& x, double& t) const;
  bool ComputeGainLoss(const Particle ptype, const size_t w,
                       std::vector<Point>& path, int& status,
                       std::vector<Seed>& secondaries,
                       const bool semiconductor = false) const;
  bool ComputeAlphaEta(const Particle ptype, std::vector<Point>& path,
                       std::vector<double>& alphas,
                       std::vector<double>& etas) const;
  bool Equilibrate(std::vector<double>& alphas) const;
  void ComputeSignal(const double q, const std::vector<Point>& path) const;
  void ComputeInducedCharge(const double q,
                            const std::vector<Point>& path) const;
  void PrintError(const std::string& fcn, const std::string& par,
                  const Particle particle,
                  const std::array<double, 3>& x) const;
};
}  // namespace Garfield
 
#endif