AvalancheMicroscopic.hh

Go to the documentation of this file.

#ifndef G_AVALANCHE_MICROSCOPIC_H
#define G_AVALANCHE_MICROSCOPIC_H
 
#include <array>
#include <cstddef>
#include <string>
#include <utility>
#include <vector>
 
#include "Garfield/GarfieldConstants.hh"
#include "Garfield/MultiProcessInterface.hh"
 
class TH1;
 
namespace Garfield {
 
class AvalancheMicroscopicGPU;
class Medium;
class ViewDrift;
class Sensor;

 
class AvalancheMicroscopic {
 public:
  AvalancheMicroscopic() : AvalancheMicroscopic(nullptr) {}
  AvalancheMicroscopic(Sensor* sensor);
  ~AvalancheMicroscopic() {}

  void SetSensor(Sensor* sensor);

  void EnablePlotting(ViewDrift* view, const std::size_t nColl = 100);
  void DisablePlotting() { m_viewer = nullptr; }
  void EnableExcitationMarkers(const bool on = true) { m_plotExcitations = on; }
  void EnableIonisationMarkers(const bool on = true) { m_plotIonisations = on; }
  void EnableAttachmentMarkers(const bool on = true) { m_plotAttachments = on; }

  void EnableSignalCalculation(const bool on = true) { m_doSignal = on; }
  void UseWeightingPotential(const bool on = true) {
    m_useWeightingPotential = on;
  }

  void EnableWeightingFieldIntegration(const bool on = true) {
    m_integrateWeightingField = on;
  }

  void UseInducedCharge(const bool on = true) { m_doInducedCharge = on; }

  void EnablePathLengthComputation(const bool on = true) {
    m_computePathLength = on;
  }

  void EnableElectronEnergyHistogramming(TH1* histo);
  void DisableElectronEnergyHistogramming() { m_histElectronEnergy = nullptr; }
  void EnableHoleEnergyHistogramming(TH1* histo);
  void DisableHoleEnergyHistogramming() { m_histHoleEnergy = nullptr; }

  void SetDistanceHistogram(TH1* histo, const char opt = 'r');
  void EnableDistanceHistogramming(const int type);
  void DisableDistanceHistogramming(const int type);
  void DisableDistanceHistogramming();
  void EnableSecondaryEnergyHistogramming(TH1* histo);
  void DisableSecondaryEnergyHistogramming() { m_histSecondary = nullptr; }

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

  void EnablePhotonTransport(const bool on = true) { m_usePhotons = on; }

  void EnableBandStructure(const bool on = true) { m_useBandStructure = on; }

  void EnableNullCollisionSteps(const bool on = true, const int nSteps = 1) {
    m_useNullCollisionSteps = on;
    m_nullCollScale = nSteps <= 1 ? 1. : 1. / nSteps;
  }

  void EnableRKNSteps(const bool on = true) { m_rknSteps = on; }
  void SetRKNTolerance(const double sTol = 1.e-10,
                       const double sMinStep = 1.e-5) {
    m_rknsteperrortol = sTol;
    m_rknMinh = sMinStep;
  }

  void SetElectronTransportCut(const double cut) { m_deltaCut = cut; }
  double GetElectronTransportCut() const { return m_deltaCut; }

  void SetPhotonTransportCut(const double cut) { m_gammaCut = cut; }
  double GetPhotonTransportCut() const { return m_gammaCut; }

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

  void EnableMagneticField(const bool on = true) {
    m_useBfieldAuto = false;
    m_useBfield = on;
  }

  void SetCollisionSteps(const unsigned int n) { m_nCollSkip = n; }

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

  void GetAvalancheSize(int& ne, int& ni) const {
    ne = m_nElectrons;
    ni = m_nIons;
  }

  std::pair<int, int> GetAvalancheSize() const {
    return std::make_pair(m_nElectrons, m_nIons);
  }
  void GetAvalancheSize(int& ne, int& nh, int& ni) const {
    ne = m_nElectrons;
    nh = m_nHoles;
    ni = m_nIons;
  }

  void GetAvalancheSizeGPU(int& ne, int& ni) const {
    ne = m_nElectronsGPU;
    ni = m_nIonsGPU;
  }
 
  struct Point {
    double x, y, z;     
    double t;           
    double energy;      
    double kx, ky, kz;  
    int band;           
  };
 
  struct Electron {
    int status = 0;           
    std::vector<Point> path;  
    size_t weight = 1;        
    double pathLength = 0.;   
  };
 
  struct Seed {
    Point pt;      
    Particle type; 
    size_t w = 1;  
  };
 
  const std::vector<Electron>& GetElectrons() const { return m_electrons; }
  const std::vector<Electron>& GetHoles() const { return m_holes; }
  size_t GetNumberOfElectronEndpoints() const { return m_electrons.size(); }
  unsigned int GetNumberOfElectronEndpointsGPU() const {
    return m_electrons_gpu.size();
  }

  void GetElectronEndpoint(const size_t i, double& x0, double& y0, double& z0,
                           double& t0, double& e0, double& x1, double& y1,
                           double& z1, double& t1, double& e1,
                           int& status) const;
  void GetElectronEndpointGPU(const size_t i, double& x0, double& y0,
                              double& z0, double& t0, double& e0, double& x1,
                              double& y1, double& z1, double& t1, double& e1,
                              int& status) const;
  size_t GetNumberOfElectronDriftLinePoints(const size_t i = 0) const;
  void GetElectronDriftLinePoint(double& x, double& y, double& z, double& t,
                                 const size_t ip, const size_t ie = 0) const;
 
  size_t GetNumberOfPhotons() const { return m_photons.size(); }
  // Status codes:
  //   -2: photon absorbed by gas molecule
  void GetPhoton(const size_t i, double& e, double& x0, double& y0, double& z0,
                 double& t0, double& x1, double& y1, double& z1, double& t1,
                 int& status) const;

  bool DriftElectron(const double x, const double y, const double z,
                     const double t, const double e, const double dx = 0.,
                     const double dy = 0., const double dz = 0.,
                     const size_t w = 1);

  bool AvalancheElectron(const double x, const double y, const double z,
                         const double t, const double e, const double dx = 0.,
                         const double dy = 0., const double dz = 0.,
                         const size_t w = 1);
  void AddElectron(const double x, const double y, const double z,
                   const double t, const double e, const double dx = 0.,
                   const double dy = 0., const double dz = 0.,
                   const size_t w = 1);
  bool ResumeAvalanche();

  void SetUserHandleStep(void (*f)(double x, double y, double z, double t,
                                   double e, double dx, double dy, double dz,
                                   bool hole));
  void UnsetUserHandleStep() { m_userHandleStep = nullptr; }
  void SetUserHandleCollision(void (*f)(double x, double y, double z, double t,
                                        int type, int level, Medium* m,
                                        double e0, double e1, double dx0,
                                        double dy0, double dz0, double dx1,
                                        double dy1, double dz1));
  void UnsetUserHandleCollision() { m_userHandleCollision = nullptr; }
  void SetUserHandleAttachment(void (*f)(double x, double y, double z, double t,
                                         int type, int level, Medium* m));
  void UnsetUserHandleAttachment() { m_userHandleAttachment = nullptr; }
  void SetUserHandleInelastic(void (*f)(double x, double y, double z, double t,
                                        int type, int level, Medium* m));
  void UnsetUserHandleInelastic() { m_userHandleInelastic = nullptr; }
  void SetUserHandleIonisation(void (*f)(double x, double y, double z, double t,
                                         int type, int level, Medium* m));
  void UnsetUserHandleIonisation() { m_userHandleIonisation = nullptr; }

  void EnableDebugging() { m_debug = true; }
  void DisableDebugging() { m_debug = false; }
 
  // class to store CPU/GPU benchmark comparisons
  struct Statistics {
    std::vector<double> gpu_stack_process_time;
    std::vector<double> cpu_stack_process_time;
    std::vector<double> gpu_stack_transport_time;
    std::vector<double> cpu_stack_transport_time;
 
    std::vector<size_t> stack_old_size;
    std::vector<size_t> stack_new_size;
  };
 
  Statistics GetStatistics() { return m_stats; }
 
  void SetRunModeOptions(MPRunMode mode, int device = -1);
  void SetMaxNumShowerLoops(int max_loops) { m_maxNumShowerLoops = max_loops; }
  void SetShowProgress(bool show_progress) { m_showProgress = show_progress; }
  void SetDebugShowerIterationAndElectronID(int iter_num, int elec_id) {
    m_debugShowerLoopNum = iter_num;
    m_debugElectronID = elec_id;
  }
 
 private:
  std::string m_className = "AvalancheMicroscopic";
 
  MPRunMode m_runMode{MPRunMode::Normal};
  int m_cudaDevice{-1};
  int m_maxNumShowerLoops{-1};
  int m_debugElectronID{-1};
  int m_debugShowerLoopNum{-1};
  bool m_showProgress{false};
  AvalancheMicroscopicGPU* m_gpuInterface{nullptr};
 
  Statistics m_stats;
 
  Sensor* m_sensor = nullptr;
 
  std::vector<Electron> m_electrons;
  std::vector<Electron> m_electrons_gpu;
  std::vector<Electron> m_holes;
 
  std::vector<Seed> m_stackStoreCPU;
  std::vector<Electron> m_stackStoreGPU;
 
  struct Photon {
    int status;             
    double energy;          
    double x0, y0, z0, t0;  
    double x1, y1, z1, t1;  
  };
  std::vector<Photon> m_photons;

  int m_nElectrons = 0;
  int m_nHoles = 0;
  int m_nIons = 0;

  int m_nElectronsGPU = 0;
  int m_nHolesGPU = 0;
  int m_nIonsGPU = 0;
 
  ViewDrift* m_viewer = nullptr;
  bool m_plotExcitations = true;
  bool m_plotIonisations = true;
  bool m_plotAttachments = true;
 
  TH1* m_histElectronEnergy = nullptr;
  TH1* m_histHoleEnergy = nullptr;
  TH1* m_histDistance = nullptr;
  char m_distanceOption = 'r';
  std::vector<int> m_distanceHistogramType;
 
  TH1* m_histSecondary = nullptr;
 
  bool m_doSignal = true;
  bool m_useWeightingPotential = true;
  bool m_integrateWeightingField = false;
  bool m_doInducedCharge = false;
 
  bool m_computePathLength = false;
  bool m_storeDriftLines = false;
  bool m_usePhotons = false;
  bool m_useBandStructure = true;
  bool m_useNullCollisionSteps = false;
  bool m_useBfieldAuto = true;
  bool m_useBfield = false;
 
  bool m_rknSteps = false;
  double m_rknsteperrortol = 1.e-10;
  double m_rknMinh = 1.e-5;
  double m_nullCollScale = 1.;
 
  // Transport cuts
  double m_deltaCut = 0.;
  double m_gammaCut = 0.;
 
  // Max. avalanche size
  unsigned int m_sizeCut = 0;
 
  size_t m_nCollSkip = 100;
  size_t m_nCollPlot = 100;
 
  bool m_hasTimeWindow = false;
  double m_tMin = 0.;
  double m_tMax = 0.;
 
  // User procedures
  void (*m_userHandleStep)(double x, double y, double z, double t, double e,
                           double dx, double dy, double dz,
                           bool hole) = nullptr;
  void (*m_userHandleCollision)(double x, double y, double z, double t,
                                int type, int level, Medium* m, double e0,
                                double e1, double dx0, double dy0, double dz0,
                                double dx1, double dy1, double dz1) = nullptr;
  void (*m_userHandleAttachment)(double x, double y, double z, double t,
                                 int type, int level, Medium* m) = nullptr;
  void (*m_userHandleInelastic)(double x, double y, double z, double t,
                                int type, int level, Medium* m) = nullptr;
  void (*m_userHandleIonisation)(double x, double y, double z, double t,
                                 int type, int level, Medium* m) = nullptr;
 
  // Switch on/off debugging messages
  bool m_debug = false;
 
  bool TransportElectrons(std::vector<Seed>& stack, const bool aval);
  int TransportElectron(
      const Seed& seed, const bool signal, 
      std::vector<double>& ts, std::vector<std::array<double, 3> >& xs,
      std::vector<Point>& path, std::vector<Seed>& stack);
  int TransportElectronBfield(
      const Seed& seed, const bool signal,
      std::vector<double>& ts, std::vector<std::array<double, 3> >& xs,
      std::vector<Point>& path, std::vector<Seed>& stack);
  int TransportElectronSc(
      const Seed& seed, const bool signal,
      std::vector<double>& ts, std::vector<std::array<double, 3> >& xs,
      std::vector<Point>& path, std::vector<Seed>& stack);
  void TransportPhoton(const double x, const double y, const double z,
                       const double t, const double e, const size_t w,
                       std::vector<Seed>& stack);
 
  bool transportParticleStack(
      const bool aval, std::vector<Seed>& stack,
      std::vector<Seed>& newParticles, const bool signal,
      const bool useBfield, const bool sc);
  void Terminate(double x0, double y0, double z0, double t0, double& x1,
                 double& y1, double& z1, double& t1) const;
 
  void CreatePenningElectron(const double x, const double y, const double z,
                             const double t, const size_t w, 
                             const double ds, const double dt, const double ep,
                             const int level, std::vector<Seed>& stack) const;
  void PlotCollision(const int cstype, const size_t did, const double x,
                     const double y, const double z, size_t& nCollPlot) const;
  void CallUserHandles(
      const int cstype, const double x, const double y, const double z, 
      const double t, const int level, Medium* medium, 
      const double en1, const double en, 
      const double kx, const double ky, const double kz,
      const double kx1, const double ky1, const double kz1) const;
  void FillDistanceHistogram(const int cstype, const double x, const double y,
                             const double z, double& xLast, double& yLast,
                             double& zLast) const;
 
 public:
  std::vector<Seed> GetStackOld() {
    return m_stackStoreCPU;
  }
  std::vector<Electron> GetStackOldGPU() { return m_stackStoreGPU; }
 
  friend class AvalancheMicroscopicGPU;
};
}  // namespace Garfield
 
#endif