Garfield++ 5.0
A toolkit for the detailed simulation of particle detectors based on ionisation measurement in gases and semiconductors
Loading...
Searching...
No Matches
Classes | Public Member Functions | List of all members
Garfield::TrackHeed Class Reference

Generate tracks using Heed++. More...

#include <TrackHeed.hh>

+ Inheritance diagram for Garfield::TrackHeed:

Classes

struct  Cluster
 
struct  Electron
 
struct  Photon
 

Public Member Functions

 TrackHeed ()
 Default constructor.
 
 TrackHeed (Sensor *sensor)
 Constructor.
 
virtual ~TrackHeed ()
 Destructor.
 
bool NewTrack (const double x0, const double y0, const double z0, const double t0, const double dx0, const double dy0, const double dz0) override
 
bool GetCluster (double &xc, double &yc, double &zc, double &tc, int &nc, double &ec, double &extra) override
 
const std::vector< Cluster > & GetClusters () const
 
bool GetCluster (double &xc, double &yc, double &zc, double &tc, int &ne, int &ni, double &ec, double &extra)
 
bool GetCluster (double &xc, double &yc, double &zc, double &tc, int &ne, int &ni, int &np, double &ec, double &extra)
 
bool GetElectron (const unsigned int i, double &x, double &y, double &z, double &t, double &e, double &dx, double &dy, double &dz)
 
bool GetIon (const unsigned int i, double &x, double &y, double &z, double &t) const
 
bool GetPhoton (const unsigned int i, double &x, double &y, double &z, double &t, double &e, double &dx, double &dy, double &dz) const
 
double GetClusterDensity () override
 
double GetStoppingPower () override
 Get the stopping power (mean energy loss [eV] per cm).
 
double GetW () const
 Return the W value of the medium (of the last simulated track).
 
double GetFanoFactor () const
 Return the Fano factor of the medium (of the last simulated track).
 
double GetPhotoAbsorptionCrossSection (const double e) const
 Return the photoabsorption cross-section at a given energy.
 
bool Initialise (Medium *medium, const bool verbose=false)
 Compute the differential cross-section for a given medium.
 
void TransportDeltaElectron (const double x0, const double y0, const double z0, const double t0, const double e0, const double dx0, const double dy0, const double dz0, int &ne, int &ni)
 
void TransportDeltaElectron (const double x0, const double y0, const double z0, const double t0, const double e0, const double dx0, const double dy0, const double dz0, int &ne)
 
void TransportPhoton (const double x0, const double y0, const double z0, const double t0, const double e0, const double dx0, const double dy0, const double dz0, int &ne, int &ni, int &np)
 
void TransportPhoton (const double x0, const double y0, const double z0, const double t0, const double e0, const double dx0, const double dy0, const double dz0, int &ne, int &ni)
 
void TransportPhoton (const double x0, const double y0, const double z0, const double t0, const double e0, const double dx0, const double dy0, const double dz0, int &ne)
 
void EnableElectricField ()
 Take the electric field into account in the stepping algorithm.
 
void DisableElectricField ()
 Do not take the electric field into account in the stepping algorithm.
 
void EnableMagneticField ()
 Take the magnetic field into account in the stepping algorithm.
 
void DisableMagneticField ()
 Do not take the magnetic field into account in the stepping algorithm.
 
void SetSteppingLimits (const double maxStep, const double radStraight, const double stepAngleStraight, const double stepAngleCurved)
 
void GetSteppingLimits (double &maxStep, double &radStraight, double &stepAngleStraight, double &stepAngleCurved)
 
void EnableCoulombScattering (const bool on=true)
 
void EnableDeltaElectronTransport ()
 Switch simulation of delta electrons on.
 
void DisableDeltaElectronTransport ()
 Switch simulation of delta electrons off.
 
void EnablePhotonReabsorption (const bool on=true)
 Simulate (or not) the photons produced in the atomic relaxation cascade.
 
void EnablePhotoAbsorptionCrossSectionOutput (const bool on)
 Write the photoabsorption cross-sections used to a text file.
 
void SetEnergyMesh (const double e0, const double e1, const int nsteps)
 
void SetParticleUser (const double m, const double z)
 
void EnableOneStepFly (const bool on)
 
- Public Member Functions inherited from Garfield::Track
 Track ()=delete
 Default constructor.
 
 Track (const std::string &name)
 Constructor.
 
virtual ~Track ()
 Destructor.
 
virtual void SetParticle (const std::string &part)
 
void SetEnergy (const double e)
 Set the particle energy.
 
void SetBetaGamma (const double bg)
 Set the relative momentum of the particle.
 
void SetBeta (const double beta)
 Set the speed ( $\beta = v/c$) of the particle.
 
void SetGamma (const double gamma)
 Set the Lorentz factor of the particle.
 
void SetMomentum (const double p)
 Set the particle momentum.
 
void SetKineticEnergy (const double ekin)
 Set the kinetic energy of the particle.
 
double GetEnergy () const
 Return the particle energy.
 
double GetBetaGamma () const
 Return the $\beta\gamma$ of the projectile.
 
double GetBeta () const
 Return the speed ( $\beta = v/c$) of the projectile.
 
double GetGamma () const
 Return the Lorentz factor of the projectile.
 
double GetMomentum () const
 Return the particle momentum.
 
double GetKineticEnergy () const
 Return the kinetic energy of the projectile.
 
double GetCharge () const
 Get the charge of the projectile.
 
double GetMass () const
 Get the mass [eV / c2] of the projectile.
 
void SetSensor (Sensor *s)
 Set the sensor through which to transport the particle.
 
void EnablePlotting (ViewDrift *viewer)
 Switch on plotting.
 
void DisablePlotting ()
 Switch off plotting.
 
void EnableDebugging ()
 Switch on debugging messages.
 
void DisableDebugging ()
 Switch off debugging messages.
 

Additional Inherited Members

- Protected Member Functions inherited from Garfield::Track
void PlotNewTrack (const double x0, const double y0, const double z0)
 
void PlotCluster (const double x0, const double y0, const double z0)
 
- Static Protected Member Functions inherited from Garfield::Track
static std::array< double, 3 > StepBfield (const double dt, const double qoverm, const double vmag, double bx, double by, double bz, std::array< double, 3 > &dir)
 
- Protected Attributes inherited from Garfield::Track
std::string m_className = "Track"
 
double m_q = -1.
 
int m_spin = 1
 
double m_mass
 
double m_energy = 0.
 
double m_beta2 = 1.
 
bool m_isElectron = false
 
std::string m_particleName = "mu-"
 
Sensor * m_sensor = nullptr
 
bool m_isChanged = true
 
ViewDrift * m_viewer = nullptr
 
bool m_debug = false
 
size_t m_plotId = 0
 

Detailed Description

Generate tracks using Heed++.

Definition at line 35 of file TrackHeed.hh.

Constructor & Destructor Documentation

◆ TrackHeed() [1/2]

Garfield::TrackHeed::TrackHeed ( )
inline

Default constructor.

Definition at line 69 of file TrackHeed.hh.

69: TrackHeed(nullptr) {}
Garfield::TrackHeed::TrackHeed
TrackHeed()
Default constructor.
Definition TrackHeed.hh:69

Referenced by TrackHeed().

◆ TrackHeed() [2/2]

Garfield::TrackHeed::TrackHeed ( Sensor * sensor)

Constructor.

Definition at line 80 of file TrackHeed.cc.

80 : Track("Heed") {
81 m_sensor = sensor;
82 m_fieldMap.reset(new Heed::HeedFieldMap());
83}
Garfield::Track::m_sensor
Sensor * m_sensor
Definition Track.hh:114
Garfield::Track::Track
Track()=delete
Default constructor.

◆ ~TrackHeed()

Garfield::TrackHeed::~TrackHeed ( )
virtual

Destructor.

Definition at line 85 of file TrackHeed.cc.

85{}

Member Function Documentation

◆ DisableDeltaElectronTransport()

void Garfield::TrackHeed::DisableDeltaElectronTransport ( )
inline

Switch simulation of delta electrons off.

Definition at line 237 of file TrackHeed.hh.

237{ m_doDeltaTransport = false; }

◆ DisableElectricField()

void Garfield::TrackHeed::DisableElectricField ( )

Do not take the electric field into account in the stepping algorithm.

Definition at line 722 of file TrackHeed.cc.

722{ m_fieldMap->UseEfield(false); }

◆ DisableMagneticField()

void Garfield::TrackHeed::DisableMagneticField ( )

Do not take the magnetic field into account in the stepping algorithm.

Definition at line 729 of file TrackHeed.cc.

729 {
730 m_fieldMap->UseBfield(false);
731 m_useBfieldAuto = false;
732}

◆ EnableCoulombScattering()

void Garfield::TrackHeed::EnableCoulombScattering ( const bool on = true)
inline

Definition at line 231 of file TrackHeed.hh.

231 {
232 m_coulombScattering = on;
233 }

◆ EnableDeltaElectronTransport()

void Garfield::TrackHeed::EnableDeltaElectronTransport ( )
inline

Switch simulation of delta electrons on.

Definition at line 235 of file TrackHeed.hh.

235{ m_doDeltaTransport = true; }

◆ EnableElectricField()

void Garfield::TrackHeed::EnableElectricField ( )

Take the electric field into account in the stepping algorithm.

Definition at line 721 of file TrackHeed.cc.

721{ m_fieldMap->UseEfield(true); }

◆ EnableMagneticField()

void Garfield::TrackHeed::EnableMagneticField ( )

Take the magnetic field into account in the stepping algorithm.

Definition at line 724 of file TrackHeed.cc.

724 {
725 m_fieldMap->UseBfield(true);
726 m_useBfieldAuto = false;
727}

Referenced by main().

◆ EnableOneStepFly()

void Garfield::TrackHeed::EnableOneStepFly ( const bool on)
inline

Definition at line 258 of file TrackHeed.hh.

258{ m_oneStepFly = on; }

◆ EnablePhotoAbsorptionCrossSectionOutput()

void Garfield::TrackHeed::EnablePhotoAbsorptionCrossSectionOutput ( const bool on)
inline

Write the photoabsorption cross-sections used to a text file.

Definition at line 245 of file TrackHeed.hh.

245 {
246 m_usePacsOutput = on;
247 }

◆ EnablePhotonReabsorption()

void Garfield::TrackHeed::EnablePhotonReabsorption ( const bool on = true)
inline

Simulate (or not) the photons produced in the atomic relaxation cascade.

Definition at line 240 of file TrackHeed.hh.

240 {
241 m_doPhotonReabsorption = on;
242 }

◆ GetCluster() [1/3]

bool Garfield::TrackHeed::GetCluster ( double & xc,
double & yc,
double & zc,
double & tc,
int & nc,
double & ec,
double & extra )
overridevirtual

Get the next "cluster" (ionising collision of the charged particle).

Parameters
xc,yc,zccoordinates of the collision
tctime of the collision
ncnumber of electrons produced
ecdeposited energy
extraadditional information (not always implemented)

Implements Garfield::Track.

Definition at line 368 of file TrackHeed.cc.

370 {
371 int ni = 0, np = 0;
372 return GetCluster(xc, yc, zc, tc, ne, ni, np, ec, extra);
373}
Garfield::TrackHeed::GetCluster
bool GetCluster(double &xc, double &yc, double &zc, double &tc, int &nc, double &ec, double &extra) override
Definition TrackHeed.cc:368

Referenced by GetCluster(), and GetCluster().

◆ GetCluster() [2/3]

bool Garfield::TrackHeed::GetCluster ( double & xc,
double & yc,
double & zc,
double & tc,
int & ne,
int & ni,
double & ec,
double & extra )

Definition at line 375 of file TrackHeed.cc.

377 {
378 int np = 0;
379 return GetCluster(xc, yc, zc, tc, ne, ni, np, ec, extra);
380}

◆ GetCluster() [3/3]

bool Garfield::TrackHeed::GetCluster ( double & xc,
double & yc,
double & zc,
double & tc,
int & ne,
int & ni,
int & np,
double & ec,
double & extra )

Get the next "cluster" (ionising collision of the charged particle).

Parameters
xc,yc,zccoordinates of the collision
tctime of the collision
nenumber of electrons
ninumber of ions
npnumber of fluorescence photons
ecdeposited energy
extraadditional information (not always implemented)

Definition at line 382 of file TrackHeed.cc.

384 {
385 // Initialise.
386 xc = yc = zc = tc = ec = extra = 0.;
387 ne = ni = np = 0;
388
389 if (m_clusters.empty()) return false;
390 // Increment the cluster index.
391 if (m_cluster < m_clusters.size()) {
392 ++m_cluster;
393 } else if (m_cluster > m_clusters.size()) {
394 m_cluster = 0;
395 }
396 if (m_cluster >= m_clusters.size()) return false;
397
398 ne = m_clusters[m_cluster].electrons.size();
399 ni = m_clusters[m_cluster].ions.size();
400 np = m_clusters[m_cluster].photons.size();
401 xc = m_clusters[m_cluster].x;
402 yc = m_clusters[m_cluster].y;
403 zc = m_clusters[m_cluster].z;
404 tc = m_clusters[m_cluster].t;
405 ec = m_clusters[m_cluster].energy;
406 extra = m_clusters[m_cluster].extra;
407 return true;
408}

◆ GetClusterDensity()

double Garfield::TrackHeed::GetClusterDensity ( )
overridevirtual

Get the cluster density (number of ionizing collisions per cm or inverse mean free path for ionization).

Reimplemented from Garfield::Track.

Definition at line 232 of file TrackHeed.cc.

232 {
233
234 if (!m_transferCs) {
235 std::cerr << m_className << "::GetClusterDensity:\n"
236 << " Ionisation cross-section is not available.\n";
237 return 0.;
238 }
239 return m_transferCs->quanC;
240}
Garfield::Track::m_className
std::string m_className
Definition Track.hh:104

◆ GetClusters()

const std::vector< Cluster > & Garfield::TrackHeed::GetClusters ( ) const
inline

Definition at line 80 of file TrackHeed.hh.

80 {
81 return m_clusters;
82 }

Referenced by main().

◆ GetElectron()

bool Garfield::TrackHeed::GetElectron ( const unsigned int i,
double & x,
double & y,
double & z,
double & t,
double & e,
double & dx,
double & dy,
double & dz )

Retrieve the properties of a conduction or delta electron in the current cluster.

Parameters
iindex of the electron
x,y,zcoordinates of the electron
ttime
ekinetic energy (only meaningful for delta-electrons)
dx,dy,dzdirection vector (only meaningful for delta-electrons)

Definition at line 410 of file TrackHeed.cc.

412 {
413
414 if (m_clusters.empty() || m_cluster >= m_clusters.size()) return false;
415 // Make sure an electron with this number exists.
416 if (i >= m_clusters[m_cluster].electrons.size()) {
417 std::cerr << m_className << "::GetElectron: Index out of range.\n";
418 return false;
419 }
420 const auto& electron = m_clusters[m_cluster].electrons[i];
421 x = electron.x;
422 y = electron.y;
423 z = electron.z;
424 t = electron.t;
425 e = electron.e;
426 dx = electron.dx;
427 dy = electron.dy;
428 dz = electron.dz;
429 return true;
430}

◆ GetFanoFactor()

double Garfield::TrackHeed::GetFanoFactor ( ) const

Return the Fano factor of the medium (of the last simulated track).

Definition at line 1124 of file TrackHeed.cc.

1124{ return m_matter->F; }

◆ GetIon()

bool Garfield::TrackHeed::GetIon ( const unsigned int i,
double & x,
double & y,
double & z,
double & t ) const

Retrieve the properties of an ion in the current cluster.

Parameters
iindex of the ion
x,y,zcoordinates of the ion
ttime

Definition at line 432 of file TrackHeed.cc.

433 {
434 if (m_clusters.empty() || m_cluster >= m_clusters.size()) return false;
435 // Make sure an ion with this index exists.
436 if (i >= m_clusters[m_cluster].ions.size()) {
437 std::cerr << m_className << "::GetIon: Index out of range.\n";
438 return false;
439 }
440 const auto& ion = m_clusters[m_cluster].ions[i];
441 x = ion.x;
442 y = ion.y;
443 z = ion.z;
444 t = ion.t;
445 return true;
446}

◆ GetPhotoAbsorptionCrossSection()

double Garfield::TrackHeed::GetPhotoAbsorptionCrossSection ( const double e) const

Return the photoabsorption cross-section at a given energy.

Definition at line 1126 of file TrackHeed.cc.

1126 {
1127
1128 if (!m_matter) return 0.;
1129 // Convert eV to MeV.
1130 const double e = 1.e-6 * en;
1131 double cs = 0.;
1132 const auto n = m_matter->apacs.size();
1133 for (size_t i = 0; i < n; ++i) {
1134 const double w = m_matter->matter->weight_quan(i);
1135 cs += m_matter->apacs[i]->get_ACS(e) * w;
1136 }
1137 // Convert Mbarn to cm-2.
1138 return cs * 1.e-18;
1139}

◆ GetPhoton()

bool Garfield::TrackHeed::GetPhoton ( const unsigned int i,
double & x,
double & y,
double & z,
double & t,
double & e,
double & dx,
double & dy,
double & dz ) const

Retrieve the properties of an unabsorbed photon.

Parameters
iindex of the photon
x,y,zcoordinates of the photon
ttime
ephoton energy
dx,dy,dzdirection vector

Definition at line 448 of file TrackHeed.cc.

450 {
451 if (m_clusters.empty() || m_cluster >= m_clusters.size()) return false;
452 // Make sure a photon with this index exists.
453 if (i >= m_clusters[m_cluster].photons.size()) {
454 std::cerr << m_className << "::GetPhoton: Index out of range.\n";
455 return false;
456 }
457 const auto& photon = m_clusters[m_cluster].photons[i];
458 x = photon.x;
459 y = photon.y;
460 z = photon.z;
461 t = photon.t;
462 e = photon.e;
463 dx = photon.dx;
464 dy = photon.dy;
465 dz = photon.dz;
466 return true;
467}

◆ GetSteppingLimits()

void Garfield::TrackHeed::GetSteppingLimits ( double & maxStep,
double & radStraight,
double & stepAngleStraight,
double & stepAngleCurved )
inline

Definition at line 223 of file TrackHeed.hh.

224 {
225 maxStep = m_maxStep;
226 radStraight = m_radStraight;
227 stepAngleStraight = m_stepAngleStraight;
228 stepAngleCurved = m_stepAngleCurved;
229 }

◆ GetStoppingPower()

double Garfield::TrackHeed::GetStoppingPower ( )
overridevirtual

Get the stopping power (mean energy loss [eV] per cm).

Reimplemented from Garfield::Track.

Definition at line 242 of file TrackHeed.cc.

242 {
243
244 if (!m_transferCs) {
245 std::cerr << m_className << "::GetStoppingPower:\n"
246 << " Ionisation cross-section is not available.\n";
247 return 0.;
248 }
249 return m_transferCs->meanC1 * 1.e6;
250}

◆ GetW()

double Garfield::TrackHeed::GetW ( ) const

Return the W value of the medium (of the last simulated track).

Definition at line 1123 of file TrackHeed.cc.

1123{ return m_matter->W * 1.e6; }

◆ Initialise()

bool Garfield::TrackHeed::Initialise ( Medium * medium,
const bool verbose = false )

Compute the differential cross-section for a given medium.

Definition at line 771 of file TrackHeed.cc.

771 {
772 // Make sure the path to the Heed database is known.
773 std::string databasePath;
774 char* dbPath = std::getenv("HEED_DATABASE");
775 if (dbPath) {
776 databasePath = dbPath;
777 } else {
778 // Try GARFIELD_INSTALL.
779 dbPath = std::getenv("GARFIELD_INSTALL");
780 if (dbPath) {
781 databasePath = std::string(dbPath) + "/share/Heed/database";
782 } else {
783 // Try GARFIELD_HOME.
784 dbPath = std::getenv("GARFIELD_HOME");
785 if (dbPath) {
786 databasePath = std::string(dbPath) + "/Heed/heed++/database";
787 } else {
788 }
789 }
790 }
791 if (databasePath.empty()) {
792 std::cerr << m_className << "::Initialise:\n"
793 << " Cannot retrieve database path (none of the"
794 << " environment variables HEED_DATABASE, GARFIELD_INSTALL,"
795 << " GARFIELD_HOME is defined).\n"
796 << " Cannot proceed.\n";
797 return false;
798 }
799 if (databasePath.back() != '/') databasePath.append("/");
800
801 if (m_debug || verbose) {
802 std::cout << m_className << "::Initialise:\n"
803 << " Database path: " << databasePath << "\n";
804 }
805
806 // Make sure the medium exists.
807 if (!medium) {
808 std::cerr << m_className << "::Initialise: Null pointer.\n";
809 return false;
810 }
811
812 // Setup the energy mesh.
813 m_energyMesh.reset(new Heed::EnergyMesh(m_emin, m_emax, m_nEnergyIntervals));
814
815 if (medium->IsGas()) {
816 if (!SetupGas(medium)) return false;
817 } else {
818 if (!SetupMaterial(medium)) return false;
819 }
820
821 // Energy transfer cross-section
822 // Set a flag indicating whether the primary particle is an electron.
823 m_transferCs.reset(new Heed::EnTransfCS(1.e-6 * m_mass, GetGamma() - 1.,
824 m_isElectron, m_matter.get(),
825 m_q));
826 if (!m_transferCs->m_ok) {
827 std::cerr << m_className << "::Initialise:\n"
828 << " Problems occured when calculating the differential"
829 << " cross-section table.\n"
830 << " Results will be inaccurate.\n";
831 }
832 if (!SetupDelta(databasePath)) return false;
833
834 if (m_debug || verbose) {
835 const double nc = m_transferCs->quanC;
836 const double dedx = m_transferCs->meanC * 1.e3;
837 const double dedx1 = m_transferCs->meanC1 * 1.e3;
838 const double w = m_matter->W * 1.e6;
839 const double f = m_matter->F;
840 const double minI = m_matter->min_ioniz_pot * 1.e6;
841 std::cout << m_className << "::Initialise:\n";
842 std::cout << " Cluster density: " << nc << " cm-1\n";
843 std::cout << " Stopping power (restricted): " << dedx << " keV/cm\n";
844 std::cout << " Stopping power (incl. tail): " << dedx1 << " keV/cm\n";
845 std::cout << " W value: " << w << " eV\n";
846 std::cout << " Fano factor: " << f << "\n";
847 std::cout << " Min. ionization potential: " << minI << " eV\n";
848 }
849
850 Heed::fixsyscoor primSys(Heed::point(0., 0., 0.), Heed::basis("primary"),
851 "primary");
852 m_chamber.reset(new HeedChamber(primSys, m_lX, m_lY, m_lZ,
853 *m_transferCs.get(), *m_deltaCs.get()));
854 m_fieldMap->SetSensor(m_sensor);
855 return true;
856}
Garfield::Track::GetGamma
double GetGamma() const
Return the Lorentz factor of the projectile.
Definition Track.hh:58
Garfield::Track::m_isElectron
bool m_isElectron
Definition Track.hh:111
Garfield::Track::m_mass
double m_mass
Definition Track.hh:108

Referenced by NewTrack(), TransportDeltaElectron(), and TransportPhoton().

◆ NewTrack()

bool Garfield::TrackHeed::NewTrack ( const double x0,
const double y0,
const double z0,
const double t0,
const double dx0,
const double dy0,
const double dz0 )
overridevirtual

Calculate a new track starting from (x0, y0, z0) at time t0 in direction (dx0, dy0, dz0).

Implements Garfield::Track.

Definition at line 87 of file TrackHeed.cc.

89 {
90 m_hasActiveTrack = false;
91
92 // Make sure the sensor has been set.
93 if (!m_sensor) {
94 std::cerr << m_className << "::NewTrack: Sensor is not defined.\n";
95 return false;
96 }
97
98 bool update = false;
99 if (!UpdateBoundingBox(update)) return false;
100
101 // Make sure the initial position is inside an ionisable medium.
102 Medium* medium = m_sensor->GetMedium(x0, y0, z0);
103 if (!medium || !medium->IsIonisable()) {
104 std::cerr << m_className << "::NewTrack:\n"
105 << " No ionisable medium at initial position.\n";
106 return false;
107 }
108
109 // Check if the medium has changed since the last call.
110 if (medium->GetName() != m_mediumName ||
111 fabs(medium->GetMassDensity() - m_mediumDensity) > 1.e-9) {
112 m_isChanged = true;
113 }
114
115 // If medium, particle or bounding box have changed,
116 // update the cross-sections.
117 if (m_isChanged) {
118 if (!Initialise(medium)) return false;
119 m_isChanged = false;
120 m_mediumName = medium->GetName();
121 m_mediumDensity = medium->GetMassDensity();
122 }
123
124 // Reset the list of clusters.
125 m_clusters.clear();
126 m_cluster = 0;
127
128 // Set the velocity vector.
129 Heed::vec velocity = NormaliseDirection(dx0, dy0, dz0);
130 velocity = velocity * Heed::CLHEP::c_light * GetBeta();
131
132 if (m_debug) {
133 std::cout << m_className << "::NewTrack:\n Track starts at (" << x0
134 << ", " << y0 << ", " << z0 << ") at time " << t0 << "\n";
135 }
136
137 // Initial position (shift with respect to bounding box center and
138 // convert from cm to mm).
139 Heed::point p0((x0 - m_cX) * 10., (y0 - m_cY) * 10., (z0 - m_cZ) * 10.);
140
141 // Setup the particle.
142 Heed::particle_def* particleType = &Heed::muon_minus_def;
143 if (m_particleName == "e-") {
144 particleType = &Heed::electron_def;
145 } else if (m_particleName == "e+") {
146 particleType = &Heed::positron_def;
147 } else if (m_particleName == "mu-") {
148 particleType = &Heed::muon_minus_def;
149 } else if (m_particleName == "mu+") {
150 particleType = &Heed::muon_plus_def;
151 } else if (m_particleName == "pi-") {
152 particleType = &Heed::pi_minus_meson_def;
153 } else if (m_particleName == "pi+") {
154 particleType = &Heed::pi_plus_meson_def;
155 } else if (m_particleName == "K-") {
156 particleType = &Heed::K_minus_meson_def;
157 } else if (m_particleName == "K+") {
158 particleType = &Heed::K_plus_meson_def;
159 } else if (m_particleName == "p") {
160 particleType = &Heed::proton_def;
161 } else if (m_particleName == "pbar") {
162 particleType = &Heed::anti_proton_def;
163 } else if (m_particleName == "d") {
164 particleType = &Heed::deuteron_def;
165 } else if (m_particleName == "alpha") {
166 particleType = &Heed::alpha_particle_def;
167 } else if (m_particleName == "exotic") {
168 // User defined particle
169 m_particle_def.reset(new Heed::particle_def(Heed::pi_plus_meson_def));
170 m_particle_def->set_mass(m_mass * 1.e-6);
171 m_particle_def->set_charge(m_q);
172 particleType = m_particle_def.get();
173 } else {
174 // Not a predefined particle, use muon definition.
175 if (m_q > 0.) {
176 particleType = &Heed::muon_minus_def;
177 } else {
178 particleType = &Heed::muon_plus_def;
179 }
180 }
181
182 Heed::HeedParticle particle(m_chamber.get(), p0, velocity, t0, particleType,
183 m_fieldMap.get(), m_coulombScattering);
184 if (m_useBfieldAuto) {
185 m_fieldMap->UseBfield(m_sensor->HasMagneticField());
186 }
187 // Set the step limits.
188 particle.set_step_limits(m_maxStep * Heed::CLHEP::cm,
189 m_radStraight * Heed::CLHEP::cm,
190 m_stepAngleStraight * Heed::CLHEP::rad,
191 m_stepAngleCurved * Heed::CLHEP::rad);
192 // Transport the particle.
193 std::vector<Heed::gparticle*> particleBank;
194 if (m_oneStepFly) {
195 particle.fly(particleBank, true);
196 } else {
197 particle.fly(particleBank);
198 }
199
200 // Sort the clusters by time.
201 std::sort(particleBank.begin(), particleBank.end(),
202 [](Heed::gparticle* p1, Heed::gparticle* p2) {
203 return p1->time() < p2->time(); });
204 // Loop over the clusters (virtual photons) created by the particle.
205 for (auto gp : particleBank) {
206 // Convert the particle to a (virtual) photon.
207 Heed::HeedPhoton* virtualPhoton = dynamic_cast<Heed::HeedPhoton*>(gp);
208 if (!virtualPhoton) {
209 std::cerr << m_className << "::NewTrack:\n"
210 << " Particle is not a virtual photon. Program bug!\n";
211 // Skip this one.
212 continue;
213 }
214 if (!AddCluster(virtualPhoton, m_clusters)) break;
215 }
216 ClearBank(particleBank);
217 Heed::gparticle::reset_counter();
218 m_cluster = m_clusters.size() + 2;
219
220 m_hasActiveTrack = true;
221
222 // Plot the track, if requested.
223 if (m_viewer) {
224 PlotNewTrack(x0, y0, z0);
225 for (const auto& cluster : m_clusters) {
226 PlotCluster(cluster.x, cluster.y, cluster.z);
227 }
228 }
229 return true;
230}
Garfield::TrackHeed::Initialise
bool Initialise(Medium *medium, const bool verbose=false)
Compute the differential cross-section for a given medium.
Definition TrackHeed.cc:771
Garfield::Track::GetBeta
double GetBeta() const
Return the speed ( ) of the projectile.
Definition Track.hh:56
Garfield::Track::m_isChanged
bool m_isChanged
Definition Track.hh:116
Garfield::Track::m_viewer
ViewDrift * m_viewer
Definition Track.hh:118
Garfield::Track::m_particleName
std::string m_particleName
Definition Track.hh:112
Garfield::Track::PlotCluster
void PlotCluster(const double x0, const double y0, const double z0)
Definition Track.cc:195
Garfield::Track::PlotNewTrack
void PlotNewTrack(const double x0, const double y0, const double z0)
Definition Track.cc:190
Heed::gparticle::reset_counter
static void reset_counter()
Reset the counter.
Definition gparticle.h:198
Heed::K_plus_meson_def
particle_def K_plus_meson_def("K_plus", "K+", 493.677 *MeV/c_squared, 1, 0.0)
Definition particle_def.h:60
Heed::alpha_particle_def
particle_def alpha_particle_def("alpha_particle", "alpha", 3727.417 *MeV/c_squared, 2 *eplus, 0.)
Definition particle_def.h:64
Heed::anti_proton_def
particle_def anti_proton_def("", "p-", proton_def)
Definition particle_def.h:55
Heed::deuteron_def
particle_def deuteron_def("deuteron", "d", 1875.613 *MeV/c_squared, eplus, 0.0)
Definition particle_def.h:63
Heed::electron_def
particle_def electron_def("electron", "e-", electron_mass_c2/c_squared, electron_charge, 0.5)
Definition particle_def.h:50
Heed::muon_minus_def
particle_def muon_minus_def("muon_minus", "mu-", 105.658367 *MeV/c_squared, electron_charge, 0.5)
Definition particle_def.h:52
Heed::fabs
DoubleAc fabs(const DoubleAc &f)
Definition DoubleAc.h:615
Heed::pi_plus_meson_def
particle_def pi_plus_meson_def("pi_plus", "pi+", 139.56755 *MeV/c_squared, eplus, 0.0)
Definition particle_def.h:58
Heed::muon_plus_def
particle_def muon_plus_def("muon_plus", "mu+", muon_minus_def)
Definition particle_def.h:53
Heed::positron_def
particle_def positron_def("positron", "e+", electron_def)
Definition particle_def.h:51
Heed::K_minus_meson_def
particle_def K_minus_meson_def("K_minus", "K-", K_plus_meson_def)
Definition particle_def.h:61
Heed::pi_minus_meson_def
particle_def pi_minus_meson_def("pi_minus", "pi-", 139.56755 *MeV/c_squared, -eplus, 0.0)
Definition particle_def.h:59
Heed::proton_def
particle_def proton_def("proton", "p+", proton_mass_c2/c_squared, eplus, 0.5)
Definition particle_def.h:54

Referenced by main().

◆ SetEnergyMesh()

void Garfield::TrackHeed::SetEnergyMesh ( const double e0,
const double e1,
const int nsteps )

Specify the energy mesh to be used.

Parameters
e0,e1lower/higher limit of the energy range [eV]
nstepsnumber of intervals

Definition at line 734 of file TrackHeed.cc.

735 {
736 if (fabs(e1 - e0) < Small) {
737 std::cerr << m_className << "::SetEnergyMesh:\n"
738 << " Invalid energy range:\n"
739 << " " << e0 << " < E [eV] < " << e1 << "\n";
740 return;
741 }
742
743 if (nsteps <= 0) {
744 std::cerr << m_className << "::SetEnergyMesh:\n"
745 << " Number of intervals must be > 0.\n";
746 return;
747 }
748
749 m_emin = 1.e-6 * std::min(e0, e1);
750 m_emax = 1.e-6 * std::max(e0, e1);
751 m_nEnergyIntervals = nsteps;
752}

◆ SetParticleUser()

void Garfield::TrackHeed::SetParticleUser ( const double m,
const double z )

Define particle mass and charge (for exotic particles). For standard particles Track::SetParticle should be used.

Definition at line 754 of file TrackHeed.cc.

754 {
755 if (fabs(z) < Small) {
756 std::cerr << m_className << "::SetParticleUser:\n"
757 << " Particle cannot have zero charge.\n";
758 return;
759 }
760 if (m < Small) {
761 std::cerr << m_className << "::SetParticleUser:\n"
762 << " Particle mass must be greater than zero.\n";
763 }
764 m_q = z;
765 m_mass = m;
766 m_isElectron = false;
767 m_spin = 0;
768 m_particleName = "exotic";
769}

◆ SetSteppingLimits()

void Garfield::TrackHeed::SetSteppingLimits ( const double maxStep,
const double radStraight,
const double stepAngleStraight,
const double stepAngleCurved )
inline

Set parameters for calculating the particle trajectory.

Parameters
maxStepmaximum step length
radStraightradius beyond which to approximate circles by polylines.
stepAngleStraightmax. angular step (in radian) when using polyline steps.
stepAngleCurvedmax. angular step (in radian) when using circular steps.

Definition at line 215 of file TrackHeed.hh.

217 {
218 m_maxStep = maxStep;
219 m_radStraight = radStraight;
220 m_stepAngleStraight = stepAngleStraight;
221 m_stepAngleCurved = stepAngleCurved;
222 }

◆ TransportDeltaElectron() [1/2]

void Garfield::TrackHeed::TransportDeltaElectron ( const double x0,
const double y0,
const double z0,
const double t0,
const double e0,
const double dx0,
const double dy0,
const double dz0,
int & ne )

Simulate a delta electron.

Parameters
x0,y0,z0initial position of the delta electron
t0initial time
e0initial kinetic energy of the delta electron
dx0,dy0,dz0initial direction of the delta electron
nenumber of electrons produced by the delta electron

Definition at line 469 of file TrackHeed.cc.

473 {
474 int ni = 0;
475 return TransportDeltaElectron(x0, y0, z0, t0, e0, dx0, dy0, dz0, ne, ni);
476}
Garfield::TrackHeed::TransportDeltaElectron
void TransportDeltaElectron(const double x0, const double y0, const double z0, const double t0, const double e0, const double dx0, const double dy0, const double dz0, int &ne, int &ni)
Definition TrackHeed.cc:478

◆ TransportDeltaElectron() [2/2]

void Garfield::TrackHeed::TransportDeltaElectron ( const double x0,
const double y0,
const double z0,
const double t0,
const double e0,
const double dx0,
const double dy0,
const double dz0,
int & ne,
int & ni )

Simulate a delta electron.

Parameters
x0,y0,z0initial position of the delta electron
t0initial time
e0initial kinetic energy of the delta electron
dx0,dy0,dz0initial direction of the delta electron
ne,ninumber of electrons/ions produced by the delta electron

Definition at line 478 of file TrackHeed.cc.

482 {
483 ne = ni = 0;
484 // Check if delta electron transport was disabled.
485 if (!m_doDeltaTransport) {
486 std::cerr << m_className << "::TransportDeltaElectron:\n"
487 << " Delta electron transport has been switched off.\n";
488 return;
489 }
490
491 // Make sure the sensor has been set.
492 if (!m_sensor) {
493 std::cerr << m_className << "::TransportDeltaElectron:\n"
494 << " Sensor is not defined.\n";
495 return;
496 }
497
498 bool update = false;
499 if (!UpdateBoundingBox(update)) return;
500
501 // Make sure the initial position is inside an ionisable medium.
502 Medium* medium = m_sensor->GetMedium(x0, y0, z0);
503 if (!medium || !medium->IsIonisable()) {
504 std::cerr << m_className << "::TransportDeltaElectron:\n"
505 << " No ionisable medium at initial position.\n";
506 return;
507 }
508
509 // Check if the medium has changed since the last call.
510 if (medium->GetName() != m_mediumName ||
511 fabs(medium->GetMassDensity() - m_mediumDensity) > 1.e-9) {
512 m_isChanged = true;
513 update = true;
514 m_hasActiveTrack = false;
515 }
516
517 // If medium or bounding box have changed, update the "chamber".
518 if (update) {
519 if (!Initialise(medium)) return;
520 m_mediumName = medium->GetName();
521 m_mediumDensity = medium->GetMassDensity();
522 }
523 m_clusters.clear();
524 m_cluster = 0;
525
526 // Initial position (shift with respect to bounding box center and
527 // convert from cm to mm).
528 Heed::point p0((x0 - m_cX) * 10., (y0 - m_cY) * 10., (z0 - m_cZ) * 10.);
529
530 Cluster cluster;
531 // Make sure the kinetic energy is positive.
532 if (e0 <= 0.) {
533 // Just create a conduction electron on the spot.
534 Electron electron;
535 electron.x = x0;
536 electron.y = y0;
537 electron.z = z0;
538 electron.t = t0;
539 cluster.electrons.push_back(std::move(electron));
540 m_clusters.push_back(std::move(cluster));
541 ne = 1;
542 return;
543 }
544
545 // Calculate the speed for the given kinetic energy.
546 const double gamma = 1. + e0 / ElectronMass;
547 const double beta = sqrt(1. - 1. / (gamma * gamma));
548 const double speed = Heed::CLHEP::c_light * beta;
549 // Set the velocity vector.
550 Heed::vec velocity = NormaliseDirection(dx0, dy0, dz0) * speed;
551
552 // Transport the electron.
553 std::vector<Heed::gparticle*> secondaries;
554 Heed::HeedDeltaElectron delta(m_chamber.get(), p0, velocity, t0, 0,
555 m_fieldMap.get());
556 delta.fly(secondaries);
557 ClearBank(secondaries);
558
559 AddElectrons(delta.conduction_electrons, cluster.electrons);
560 AddElectrons(delta.conduction_ions, cluster.ions);
561 ne = cluster.electrons.size();
562 ni = cluster.ions.size();
563 m_clusters.push_back(std::move(cluster));
564}
Heed::sqrt
DoubleAc sqrt(const DoubleAc &f)
Definition DoubleAc.cpp:314

Referenced by TransportDeltaElectron().

◆ TransportPhoton() [1/3]

void Garfield::TrackHeed::TransportPhoton ( const double x0,
const double y0,
const double z0,
const double t0,
const double e0,
const double dx0,
const double dy0,
const double dz0,
int & ne )

Simulate a photon.

Parameters
x0,y0,z0initial position of the photon
t0initial time
e0initial energy of the photon
dx0,dy0,dz0initial direction of the photon
nenumber of electrons produced by the photon

Definition at line 566 of file TrackHeed.cc.

569 {
570 int ni = 0, np = 0;
571 TransportPhoton(x0, y0, z0, t0, e0, dx0, dy0, dz0, ne, ni, np);
572}
Garfield::TrackHeed::TransportPhoton
void TransportPhoton(const double x0, const double y0, const double z0, const double t0, const double e0, const double dx0, const double dy0, const double dz0, int &ne, int &ni, int &np)
Definition TrackHeed.cc:583

◆ TransportPhoton() [2/3]

void Garfield::TrackHeed::TransportPhoton ( const double x0,
const double y0,
const double z0,
const double t0,
const double e0,
const double dx0,
const double dy0,
const double dz0,
int & ne,
int & ni )

Simulate a photon.

Parameters
x0,y0,z0initial position of the photon
t0initial time
e0initial energy of the photon
dx0,dy0,dz0initial direction of the photon
nenumber of electrons produced by the photon
ninumber of ions produced by the photon

Definition at line 574 of file TrackHeed.cc.

578 {
579 int np = 0;
580 TransportPhoton(x0, y0, z0, t0, e0, dx0, dy0, dz0, ne, ni, np);
581}

◆ TransportPhoton() [3/3]

void Garfield::TrackHeed::TransportPhoton ( const double x0,
const double y0,
const double z0,
const double t0,
const double e0,
const double dx0,
const double dy0,
const double dz0,
int & ne,
int & ni,
int & np )

Simulate a photon.

Parameters
x0,y0,z0initial position of the photon
t0initial time
e0initial energy of the photon
dx0,dy0,dz0initial direction of the photon
nenumber of electrons produced by the photon
ninumber of ions produced by the photon
npnumber of fluorescence photons

Definition at line 583 of file TrackHeed.cc.

587 {
588 ne = ni = np = 0;
589 // Make sure the energy is positive.
590 if (e0 <= 0.) {
591 std::cerr << m_className << "::TransportPhoton:\n"
592 << " Photon energy must be positive.\n";
593 return;
594 }
595
596 // Make sure the sensor has been set.
597 if (!m_sensor) {
598 std::cerr << m_className << "::TransportPhoton: Sensor is not defined.\n";
599 return;
600 }
601
602 bool update = false;
603 if (!UpdateBoundingBox(update)) return;
604
605 // Make sure the initial position is inside an ionisable medium.
606 Medium* medium = m_sensor->GetMedium(x0, y0, z0);
607 if (!medium || !medium->IsIonisable()) {
608 std::cerr << m_className << "::TransportPhoton:\n"
609 << " No ionisable medium at initial position.\n";
610 return;
611 }
612
613 // Check if the medium has changed since the last call.
614 if (medium->GetName() != m_mediumName ||
615 fabs(medium->GetMassDensity() - m_mediumDensity) > 1.e-9) {
616 m_isChanged = true;
617 update = true;
618 }
619
620 // If medium or bounding box have changed, update the "chamber".
621 if (update) {
622 if (!Initialise(medium)) return;
623 m_mediumName = medium->GetName();
624 m_mediumDensity = medium->GetMassDensity();
625 }
626
627 // Clusters from the current track will be lost.
628 m_hasActiveTrack = false;
629 m_clusters.clear();
630 m_cluster = 0;
631 Cluster cluster;
632
633 // Set the direction vector.
634 Heed::vec velocity = NormaliseDirection(dx0, dy0, dz0);
635 velocity = velocity * Heed::CLHEP::c_light;
636
637 // Initial position (shift with respect to bounding box center and
638 // convert from cm to mm).
639 Heed::point p0((x0 - m_cX) * 10., (y0 - m_cY) * 10., (z0 - m_cZ) * 10.);
640
641 // Create and transport the photon.
642 Heed::HeedPhoton photon(m_chamber.get(), p0, velocity, t0, 0, e0 * 1.e-6,
643 m_fieldMap.get());
644 std::vector<Heed::gparticle*> secondaries;
645 photon.fly(secondaries);
646 if (secondaries.empty()) {
647 Photon unabsorbedPhoton;
648 unabsorbedPhoton.x = photon.position().x * 0.1 + m_cX;
649 unabsorbedPhoton.y = photon.position().y * 0.1 + m_cY;
650 unabsorbedPhoton.z = photon.position().z * 0.1 + m_cZ;
651 unabsorbedPhoton.t = photon.time();
652 unabsorbedPhoton.e = photon.m_energy * 1.e6;
653 unabsorbedPhoton.dx = photon.direction().x;
654 unabsorbedPhoton.dy = photon.direction().y;
655 unabsorbedPhoton.dz = photon.direction().z;
656 cluster.photons.push_back(std::move(unabsorbedPhoton));
657 }
658
659 while (!secondaries.empty()) {
660 std::vector<Heed::gparticle*> newSecondaries;
661 // Loop over the particle bank and look for daughter particles.
662 for (auto gp : secondaries) {
663 // Is it a delta electron?
664 auto delta = dynamic_cast<Heed::HeedDeltaElectron*>(gp);
665 if (delta) {
666 if (m_doDeltaTransport) {
667 // Transport the delta electron.
668 delta->fly(newSecondaries);
669 // Add the conduction electrons and ions to the list.
670 AddElectrons(delta->conduction_electrons, cluster.electrons);
671 AddElectrons(delta->conduction_ions, cluster.ions);
672 } else {
673 // Add the delta electron to the list, for later use.
674 Electron deltaElectron;
675 deltaElectron.x = delta->position().x * 0.1 + m_cX;
676 deltaElectron.y = delta->position().y * 0.1 + m_cY;
677 deltaElectron.z = delta->position().z * 0.1 + m_cZ;
678 deltaElectron.t = delta->time();
679 deltaElectron.e = delta->kinetic_energy() * 1.e6;
680 deltaElectron.dx = delta->direction().x;
681 deltaElectron.dy = delta->direction().y;
682 deltaElectron.dz = delta->direction().z;
683 cluster.electrons.push_back(std::move(deltaElectron));
684 }
685 continue;
686 }
687 // Check if it is a fluorescence photon.
688 auto fluorescencePhoton = dynamic_cast<Heed::HeedPhoton*>(gp);
689 if (!fluorescencePhoton) {
690 std::cerr << m_className << "::TransportPhoton:\n"
691 << " Unknown secondary particle.\n";
692 ClearBank(secondaries);
693 ClearBank(newSecondaries);
694 return;
695 }
696 if (m_doPhotonReabsorption) {
697 fluorescencePhoton->fly(newSecondaries);
698 } else {
699 Photon unabsorbedPhoton;
700 unabsorbedPhoton.x = fluorescencePhoton->position().x * 0.1 + m_cX;
701 unabsorbedPhoton.y = fluorescencePhoton->position().y * 0.1 + m_cY;
702 unabsorbedPhoton.z = fluorescencePhoton->position().z * 0.1 + m_cZ;
703 unabsorbedPhoton.t = fluorescencePhoton->time();
704 unabsorbedPhoton.e = fluorescencePhoton->m_energy * 1.e6;
705 unabsorbedPhoton.dx = fluorescencePhoton->direction().x;
706 unabsorbedPhoton.dy = fluorescencePhoton->direction().y;
707 unabsorbedPhoton.dz = fluorescencePhoton->direction().z;
708 cluster.photons.push_back(std::move(unabsorbedPhoton));
709 }
710 }
711 secondaries.swap(newSecondaries);
712 ClearBank(newSecondaries);
713 }
714 ClearBank(secondaries);
715 ne = cluster.electrons.size();
716 ni = cluster.ions.size();
717 np = cluster.photons.size();
718 m_clusters.push_back(std::move(cluster));
719}

Referenced by TransportPhoton(), and TransportPhoton().

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