AvalancheGridSpaceCharge.hh

Go to the documentation of this file.

#ifndef GARFIELD_AVALANCHEGRIDSPACECHARGE_HH
#define GARFIELD_AVALANCHEGRIDSPACECHARGE_HH
 
#include <algorithm>
#include <string>
#include <utility>
#include <vector>
#include <array>
 
namespace Garfield {
class Sensor;
class AvalancheMicroscopic;
class ComponentParallelPlate;

class AvalancheGridSpaceCharge {
 public:
  AvalancheGridSpaceCharge() : AvalancheGridSpaceCharge(nullptr) {}
  AvalancheGridSpaceCharge(Sensor *sensor);
  ~AvalancheGridSpaceCharge() = default;

  void Reset();

  void EnableDebugging(const bool option = true) { m_bDebug = option; }

  void EnableStickyAnode(const bool option = true) { m_bStick = option; }

  void EnableTOF(const bool option = true) { m_bUseTOF = option; }

  void EnableDiffusion(const bool option = true) { m_bDiffusion = option; }

  void EnableSpaceChargeEffect(const bool option = true) {
    m_bSpaceCharge = option;
  }

  void EnableAdaptiveTimeStepping(const bool option = true) {
    m_bAdaptiveTime = option;
  }

  void EnableMC(const bool option = true) { m_bMC = option; }

  void SetNCrit(const long NCrit = 1e8) { m_lNCrit = NCrit; }

  void SetFieldCalculation(const std::string &option = "coulomb",
                           const int nof_approx = 1.) {
    m_sFieldOption = std::move(option);
    m_iFieldApprox = std::move(nof_approx);
  }

  void SetK(float option = 0.95) { m_fStreamerK = option; }

  void SetStopAtK(bool option = true) { m_bStopAtK = option; }

  void SetSensor(Sensor *sensor);

  void Set2dGrid(double zmin, double zmax, int zsteps, double rmax, int rsteps);

  void AddElectrons(AvalancheMicroscopic *avmc);

  void AddElectron(double x, double y, double z, double t = 0, int n = 1);

  void AddExtraElectron(double y, int n = 1);

  void StartGridAvalanche(double dtime = -1);

  long GetAvalancheSize() const { return m_nTotPosIons; }

  double GetMeanDistance();

  [[nodiscard]] long ReachedKPercent() const {
    if (m_bFieldK)
      return m_lElectronsK;
    else
      return -1;
  }

  [[nodiscard]] const std::vector<std::pair<double, long>>
  GetElectronEvolution() const {
    return m_vNElectronEvolution;
  }

  void ExportGrid(const std::string &filename);
 
 private:
  struct GridNode {
    long nElectron = 0;  
    double nPosIon = 0;  
    double nNegIon = 0;  
    // holder memories for stepping in time:
    long nElectronHolder = 0;  
    double nPosIonHolder = 0;  
    double nNegIonHolder = 0;  
 
    double townsend = 0;    
    double attachment = 0;  
    double velocity = 0.;
    double dSigmaL = 0;
    double dSigmaT = 0;
 
    double Wv = 0;  
    double Wr = 0;  
    double townsendPT = 0;
    double attachmentPT = 0;
    double eFieldR = 0;
    double eFieldZ = 0;
 
    double time = 0.;  
 
    bool anode = false;  
    int layerIndex = 0;
    int gasGapIndex = 0;
    bool isGasGap = true;
  };
 
  struct Point {
    double x, y, z;  
    double t;        
 
    int gasLayerIndex;
  };
 
  // Prepare grid and place stored electrons from AvalancheMicroscopic import
  void PrepareElectronsFromMicroscopicAvalanche();
 
  // Assign electron to the closest grid point
  bool SnapTo2dGrid(double x, double y, double z, long n = 1, int gasLayer = 0);
 
  // Prepare the mesh with the ComponentParallelPlate
  void Prepare2dMesh();
 
  // Transports the electrons/nodes a timestep
  bool TransportTimeStep();
 
  // Diffuses the electrons/nodes a timestep
  void DiffuseTimeStep(double dx, long nElectron, double nPosIon,
                       double nNegIon, int iz, int ir, int gasGap);
 
  // Redistributes the charges
  void DistributeCharges(long nElectron, double nPosIon, double nNegIon, int iz,
                         int ir, double stepZ, double stepR, int gasGap);
 
  // Calculate the field from all the contributions to the bin of interest. May
  // need much more functionalities/tables.
  void GetLocalField(int iz, int ir, double &eFieldZ, double &eFieldR,
                     const std::string &fieldOption, int gasGap);
 
  // Calculate the field of charged ring in vacuum using coulomb potentials and
  // indices
  void GetFreeChargedRing(int iz, int ir, int fz, int fr, double &eFieldZ,
                          double &eFieldR);
 
  // Calculate the field of charged ring in vacuum using coulomb potentials and
  // coordinates
  void GetFreeChargedRing(double zi, double ri, double zf, double rf,
                          double &eFieldZ, double &eFieldR);
 
  // Get field at (zi, ri) from N charges at (zf, rf) either as a ring or a
  // coulomb ball (rf = 0) if i and f are too close it is considered as self
  // interaction and not included
  bool AddFieldFromChargeAt(int iz, int ir, int fz, int fr, double N,
                            double &eFieldZ, double &eFieldR);
 
  // Get field at (zi, ri) from N charges at (zf, rf) either as a ring or a
  // coulomb ball (rf = 0) if i and f are too close it is considered as self
  // interaction and not included
  bool AddFieldFromChargeAt(int iz, int ir, double zf, double rf, double N,
                            double &eFieldZ, double &eFieldR);
 
  // Get swarm parameters at electric field magnitude
  void GetSwarmParameters(double MagEField, double &alpha, double &eta,
                          double &drift, double &dSigmaL, double &dSigmaT,
                          double &wv, double &wr, double &alphaPT,
                          double &etaPT, int gasGap);
 
  // Change from 2dGrid to Global coordinates
  void GetGlobalCoordinates(double r, double z, double phi, double &xg,
                            double &yg, double &zg, int gasGap);
 
  // Import elliptic integral values
  void ImportEllipticIntegralValues(const std::string &filename);
 
  // Gets elliptic integrals via list
  void GetEllipticIntegrals(double x, double &K, double &E);
 
  // Get from index the gas gap number, else -1
  int GetGasGapNumber(int layerIndex) {
    auto it =
        std::find(m_vIndexGasGaps.begin(), m_vIndexGasGaps.end(), layerIndex);
    return (it != m_vIndexGasGaps.end())
               ? std::distance(m_vIndexGasGaps.begin(), it)
               : -1;
  }
 
 private:
  std::string m_className = "AvalancheGridSpaceCharge";
 
  bool m_bDebug = false;
  bool m_bDiffusion = false;
  bool m_bStick = true;
  // boolean for AvalancheElectron
  bool m_bDriftAvalanche = false;
  // boolean for ImportElectronsFromAvalancheMicroscopic
  bool m_bImportAvalanche = false;
  bool m_bPreparedImportAvalanche = false;
  long m_lNCrit = 1e8;
  bool m_bSpaceCharge = true;
 
  float m_fStreamerK = 0.95;
  bool m_bStopAtK = false;
  bool m_bFieldK = false;
  long m_lElectronsK;
 
  bool m_bAdaptiveTime = true;
  bool m_bImportElliptic = false;
  bool m_bUseTOF = true;
  bool m_bWrAvailable = true;
  bool m_bRatesAvailable = true;
 
  bool m_bMC = true;
 
  int m_iFieldApprox = 1;    //< order of approximation in Set(1,2,3,...)
  double m_dMinGroups = 50;  //< same values as lippmann
 
  ComponentParallelPlate *m_pp = nullptr;
  Sensor *m_sensor = nullptr;
 
  std::vector<double> m_zGrid;  
  int m_zSteps = 0;             
  double m_zStepSize = 0.;      
 
  std::vector<double> m_rGrid;  
  int m_rSteps = 0.;            
  double m_rStepSize = 0.;      
 
  bool m_isgridset = false;  
  long m_nTotElectron = 0;   
  long m_nTotPosIons = 0;    
 
  double m_time = 0.;   
  double m_time0 = 0.;  
  double m_dt = 0.;     
  bool m_run = true;
 
  std::vector<std::vector<int>>
      m_zGasGapBoundaries;  
 
  std::vector<std::vector<GridNode>> m_grid;  
  std::vector<std::vector<Point>> m_vElectrons;
 
  std::vector<std::pair<double, long>> m_vNElectronEvolution;
  std::vector<long> m_vGroupSizes = {
      1500, 800, 400, 200, 100,
      50,   20,  10,  5,   2};  

  std::vector<int> m_vIndexGasGaps = {0};

  std::vector<std::vector<double>> m_vCoNGasLayer{};
  std::vector<double> m_vYPointInGasGap{};

  std::vector<double> m_ezBkg = {0};
  std::vector<int> m_vSaturatedGaps{};
 
  std::string m_sFieldOption = "coulomb";
  enum class Elliptic : std::size_t
  {
    X,
    K,
    E
  };
  static const constexpr std::size_t elliptic_size{29981};
  static const std::array<std::array<double,3>, elliptic_size> m_elliptic; 
};
 
}  // namespace Garfield
 
#endif  // GARFIELD_AVALANCHEGRIDSPACECHARGE_HH