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

#include <TrackBichsel.hh>

+ Inheritance diagram for Garfield::TrackBichsel:

Public Member Functions

 TrackBichsel ()
 Constructor.
 
virtual ~TrackBichsel ()
 Destructor.
 
virtual bool NewTrack (const double x0, const double y0, const double z0, const double t0, const double dx0, const double dy0, const double dz0)
 
virtual bool GetCluster (double &xcls, double &ycls, double &zcls, double &tcls, int &n, double &e, double &extra)
 
virtual double GetClusterDensity ()
 
virtual double GetStoppingPower ()
 Get the stopping power (mean energy loss [eV] per cm).
 
void SetDataFile (const std::string &filename)
 
- Public Member Functions inherited from Garfield::Track
 Track ()
 Constructor.
 
virtual ~Track ()
 Destructor.
 
virtual void SetParticle (const std::string &part)
 Set the type of particle.
 
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
 
double GetBetaGamma () const
 
double GetBeta () const
 
double GetGamma () const
 
double GetMomentum () const
 
double GetKineticEnergy () const
 
double GetCharge () const
 Get the charge of the projectile.
 
double GetMass () const
 Get the mass [eV / c2] of the projectile.
 
void SetSensor (Sensor *s)
 
virtual bool NewTrack (const double x0, const double y0, const double z0, const double t0, const double dx0, const double dy0, const double dz0)=0
 
virtual bool GetCluster (double &xcls, double &ycls, double &zcls, double &tcls, int &n, double &e, double &extra)=0
 
virtual double GetClusterDensity ()
 
virtual double GetStoppingPower ()
 Get the stopping power (mean energy loss [eV] per cm).
 
void EnablePlotting (ViewDrift *viewer)
 
void DisablePlotting ()
 
void EnableDebugging ()
 
void DisableDebugging ()
 

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)
 
- Protected Attributes inherited from Garfield::Track
std::string m_className
 
double m_q
 
int m_spin
 
double m_mass
 
double m_energy
 
double m_beta2
 
bool m_isElectron
 
std::string m_particleName
 
Sensor * m_sensor
 
bool m_isChanged
 
bool m_usePlotting
 
ViewDrift * m_viewer
 
bool m_debug
 
int m_plotId
 

Detailed Description

Generate tracks using differential cross-sections for silicon computed by Hans Bichsel. References:

Definition at line 14 of file TrackBichsel.hh.

Constructor & Destructor Documentation

◆ TrackBichsel()

Garfield::TrackBichsel::TrackBichsel ( )

Constructor.

Definition at line 15 of file TrackBichsel.cc.

16 : m_bg(3.16228),
17 m_speed(SpeedOfLight * m_bg / sqrt(1. + m_bg * m_bg)),
18 m_x(0.),
19 m_y(0.),
20 m_z(0.),
21 m_t(0.),
22 m_dx(0.),
23 m_dy(0.),
24 m_dz(1.),
25 m_imfp(4.05090e4),
26 m_datafile("SiM0invw.inv"),
27 m_iCdf(2),
28 m_nCdfEntries(-1),
29 m_isInitialised(false),
30 m_isInMedium(false) {
31
32 m_className = "TrackBichsel";
33}
Garfield::Track::m_className
std::string m_className
Definition: Track.hh:80
Heed::sqrt
DoubleAc sqrt(const DoubleAc &f)
Definition: DoubleAc.cpp:314

◆ ~TrackBichsel()

virtual Garfield::TrackBichsel::~TrackBichsel ( )
inlinevirtual

Destructor.

Definition at line 20 of file TrackBichsel.hh.

20{}

Member Function Documentation

◆ GetCluster()

bool Garfield::TrackBichsel::GetCluster ( double &  xcls,
double &  ycls,
double &  zcls,
double &  tcls,
int &  n,
double &  e,
double &  extra 
)
virtual

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

Parameters
xcls,ycls,zclscoordinates of the collision
tclstime of the collision
nnumber of electrons produced
edeposited energy
extraadditional information (not always implemented)

Implements Garfield::Track.

Definition at line 111 of file TrackBichsel.cc.

112 {
113
114 if (!m_isInitialised || !m_isInMedium) return false;
115
116 const double d = -log(RndmUniformPos()) / m_imfp;
117 m_x += m_dx * d;
118 m_y += m_dy * d;
119 m_z += m_dz * d;
120 m_t += d / m_speed;
121
122 xcls = m_x;
123 ycls = m_y;
124 zcls = m_z;
125 tcls = m_t;
126 n = 0;
127 e = 0.;
128 extra = 0.;
129
130 Medium* medium;
131 if (!m_sensor->GetMedium(m_x, m_y, m_z, medium)) {
132 m_isInMedium = false;
133 if (m_debug) {
134 std::cout << m_className << "::GetCluster: Particle left the medium.\n";
135 }
136 return false;
137 }
138
139 if (medium->GetName() != "Si" || !medium->IsIonisable()) {
140 m_isInMedium = false;
141 if (m_debug) {
142 std::cout << m_className << "::GetCluster: Particle left the medium.\n";
143 }
144 return false;
145 }
146
147 const double u = m_nCdfEntries * RndmUniform();
148 const int j = int(u);
149 if (j == 0) {
150 e = 0. + u * m_cdf[0][m_iCdf];
151 } else if (j >= m_nCdfEntries) {
152 e = m_cdf[m_nCdfEntries - 1][m_iCdf];
153 } else {
154 e = m_cdf[j - 1][m_iCdf] +
155 (u - j) * (m_cdf[j][m_iCdf] - m_cdf[j - 1][m_iCdf]);
156 }
157
158 return true;
159}
Garfield::Sensor::GetMedium
bool GetMedium(const double x, const double y, const double z, Medium *&medium)
Get the medium at (x, y, z).
Definition: Sensor.cc:150
Garfield::Track::m_sensor
Sensor * m_sensor
Definition: Track.hh:90
Garfield::Track::m_debug
bool m_debug
Definition: Track.hh:97
Garfield::RndmUniform
double RndmUniform()
Draw a random number uniformly distributed in the range [0, 1).
Definition: Random.hh:14
Garfield::RndmUniformPos
double RndmUniformPos()
Draw a random number uniformly distributed in the range (0, 1).
Definition: Random.hh:17

◆ GetClusterDensity()

double Garfield::TrackBichsel::GetClusterDensity ( )
virtual

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

Reimplemented from Garfield::Track.

Definition at line 161 of file TrackBichsel.cc.

161 {
162
163 const int nEntries = 38;
164
165 const double tabBg[nEntries] = {
166 0.316, 0.398, 0.501, 0.631, 0.794, 1.000, 1.259, 1.585,
167 1.995, 2.512, 3.162, 3.981, 5.012, 6.310, 7.943, 10.000,
168 12.589, 15.849, 19.953, 25.119, 31.623, 39.811, 50.119, 63.096,
169 79.433, 100.000, 125.893, 158.489, 199.526, 251.189, 316.228, 398.107,
170 501.187, 630.958, 794.329, 1000.000, 1258.926, 1584.894};
171
172 const double tabImfp[nEntries] = {
173 30.32496, 21.14965, 15.06555, 11.05635, 8.43259, 6.72876, 5.63184,
174 4.93252, 4.49174, 4.21786, 4.05090, 3.95186, 3.89531, 3.86471,
175 3.84930, 3.84226, 3.83952, 3.83887, 3.83912, 3.83970, 3.84035,
176 3.84095, 3.84147, 3.84189, 3.84223, 3.84249, 3.84269, 3.84283,
177 3.84293, 3.84300, 3.84304, 3.84308, 3.84310, 3.84311, 3.84312,
178 3.84313, 3.84313, 3.84314};
179
180 if (m_isChanged) m_bg = GetBetaGamma();
181
182 if (m_bg < tabBg[0]) {
183 if (m_debug) {
184 std::cerr << m_className << "::GetClusterDensity:\n"
185 << " Bg is below the tabulated range.\n";
186 }
187 return tabImfp[0] * 1.e4;
188 } else if (m_bg > tabBg[nEntries - 1]) {
189 return tabImfp[nEntries - 1] * 1.e4;
190 }
191
192 // Locate the requested energy in the table
193 int iLow = 0;
194 int iUp = nEntries - 1;
195 while (iUp - iLow > 1) {
196 const int iM = (iUp + iLow) >> 1;
197 if (m_bg >= tabBg[iM]) {
198 iLow = iM;
199 } else {
200 iUp = iM;
201 }
202 }
203
204 if (fabs(m_bg - tabBg[iLow]) < 1.e-6 * (tabBg[iUp] - tabBg[iLow])) {
205 return tabImfp[iLow] * 1.e4;
206 }
207 if (fabs(m_bg - tabBg[iUp]) < 1.e-6 * (tabBg[iUp] - tabBg[iLow])) {
208 return tabImfp[iUp] * 1.e4;
209 }
210
211 // Log-log interpolation
212 const double logX0 = log(tabBg[iLow]);
213 const double logX1 = log(tabBg[iUp]);
214 const double logY0 = log(tabImfp[iLow]);
215 const double logY1 = log(tabImfp[iUp]);
216 double d = logY0 + (log(m_bg) - logX0) * (logY1 - logY0) / (logX1 - logX0);
217 return 1.e4 * exp(d);
218}
Garfield::Track::GetBetaGamma
double GetBetaGamma() const
Definition: Track.hh:39
Garfield::Track::m_isChanged
bool m_isChanged
Definition: Track.hh:92
Heed::exp
DoubleAc exp(const DoubleAc &f)
Definition: DoubleAc.cpp:377
Heed::fabs
DoubleAc fabs(const DoubleAc &f)
Definition: DoubleAc.h:615

Referenced by NewTrack().

◆ GetStoppingPower()

double Garfield::TrackBichsel::GetStoppingPower ( )
virtual

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

Reimplemented from Garfield::Track.

Definition at line 220 of file TrackBichsel.cc.

220 {
221
222 const unsigned int nEntries = 51;
223
224 const double tabBg[nEntries] = {
225 0.316, 0.398, 0.501, 0.631, 0.794, 1.000, 1.259,
226 1.585, 1.995, 2.512, 3.162, 3.981, 5.012, 6.310,
227 7.943, 10.000, 12.589, 15.849, 19.953, 25.119, 31.623,
228 39.811, 50.119, 63.096, 79.433, 100.000, 125.893, 158.489,
229 199.526, 251.189, 316.228, 398.107, 501.187, 630.958, 794.329,
230 1000.000, 1258.926, 1584.894, 1995.263, 2511.888, 3162.280, 3981.074,
231 5011.875, 6309.578, 7943.287, 10000.010, 12589.260, 15848.940, 19952.640,
232 25118.880, 31622.800};
233
234 const double tabdEdx[nEntries] = {
235 2443.71800, 1731.65600, 1250.93400, 928.69920, 716.37140, 578.28850,
236 490.83670, 437.33820, 406.58490, 390.95170, 385.29000, 386.12000,
237 391.07730, 398.53930, 407.39420, 416.90860, 426.63010, 436.30240,
238 445.78980, 455.02530, 463.97370, 472.61410, 480.92980, 488.90240,
239 496.51900, 503.77130, 510.65970, 517.19570, 523.39830, 529.29120,
240 534.90670, 540.27590, 545.42880, 550.39890, 555.20800, 559.88820,
241 564.45780, 568.93850, 573.34700, 577.69140, 581.99010, 586.25090,
242 590.47720, 594.68660, 598.86880, 603.03510, 607.18890, 611.33250,
243 615.46810, 619.59740, 623.72150};
244
245 if (m_isChanged) m_bg = GetBetaGamma();
246
247 if (m_bg < tabBg[0]) {
248 if (m_debug) {
249 std::cerr << m_className << "::GetStoppingPower:\n";
250 std::cerr << " Bg is below the tabulated range.\n";
251 }
252 return tabdEdx[0] * 1.e4;
253 } else if (m_bg > tabBg[nEntries - 1]) {
254 return tabdEdx[nEntries - 1] * 1.e4;
255 }
256
257 // Locate the requested energy in the table
258 int iLow = 0;
259 int iUp = nEntries - 1;
260 int iM;
261 while (iUp - iLow > 1) {
262 iM = (iUp + iLow) >> 1;
263 if (m_bg >= tabBg[iM]) {
264 iLow = iM;
265 } else {
266 iUp = iM;
267 }
268 }
269
270 if (m_debug) {
271 std::cout << m_className << "::GetStoppingPower:\n";
272 std::cout << " Bg = " << m_bg << "\n";
273 std::cout << " Interpolating between " << tabBg[iLow] << " and "
274 << tabBg[iUp] << "\n";
275 }
276
277 if (fabs(m_bg - tabBg[iLow]) < 1.e-6 * (tabBg[iUp] - tabBg[iLow])) {
278 return tabdEdx[iLow] * 1.e4;
279 }
280 if (fabs(m_bg - tabBg[iUp]) < 1.e-6 * (tabBg[iUp] - tabBg[iLow])) {
281 return tabdEdx[iUp] * 1.e4;
282 }
283
284 // Log-log interpolation
285 const double logX0 = log(tabBg[iLow]);
286 const double logX1 = log(tabBg[iUp]);
287 const double logY0 = log(tabdEdx[iLow]);
288 const double logY1 = log(tabdEdx[iUp]);
289 const double dedx =
290 logY0 + (log(m_bg) - logX0) * (logY1 - logY0) / (logX1 - logX0);
291 return 1.e4 * exp(dedx);
292}

◆ NewTrack()

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

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

Implements Garfield::Track.

Definition at line 35 of file TrackBichsel.cc.

37 {
38
39 // Make sure a sensor has been defined.
40 if (!m_sensor) {
41 std::cerr << m_className << "::NewTrack:\n"
42 << " Sensor is not defined.\n";
43 m_isInMedium = false;
44 return false;
45 }
46
47 // If not yet done, load the cross-section table from file.
48 if (!m_isInitialised) {
49 if (!LoadCrossSectionTable(m_datafile)) {
50 std::cerr << m_className << "::NewTrack:\n"
51 << " Cross-section table could not be loaded.\n";
52 return false;
53 }
54 m_isInitialised = true;
55 }
56
57 // Make sure we are inside a medium.
58 Medium* medium;
59 if (!m_sensor->GetMedium(x0, y0, z0, medium)) {
60 std::cerr << m_className << "::NewTrack:\n"
61 << " No medium at initial position.\n";
62 m_isInMedium = false;
63 return false;
64 }
65
66 // Check if the medium is silicon.
67 if (medium->GetName() != "Si") {
68 std::cerr << m_className << "::NewTrack:\n"
69 << " Medium at initial position is not silicon.\n";
70 m_isInMedium = false;
71 return false;
72 }
73
74 // Check if primary ionisation has been enabled.
75 if (!medium->IsIonisable()) {
76 std::cerr << m_className << "::NewTrack:\n"
77 << " Medium at initial position is not ionisable.\n";
78 m_isInMedium = false;
79 return false;
80 }
81
82 m_isInMedium = true;
83 m_x = x0;
84 m_y = y0;
85 m_z = z0;
86 m_t = t0;
87
88 // Normalise the direction vector.
89 const double d = sqrt(dx0 * dx0 + dy0 * dy0 + dz0 * dz0);
90 if (d < Small) {
91 // In case of a null vector, choose a random direction.
92 RndmDirection(m_dx, m_dy, m_dz);
93 } else {
94 m_dx = dx0 / d;
95 m_dy = dy0 / d;
96 m_dz = dz0 / d;
97 }
98
99 // If the particle properties have changed, update the cross-section table.
100 if (m_isChanged) {
101 m_bg = GetBetaGamma();
102 m_imfp = GetClusterDensity();
103 m_speed = SpeedOfLight * GetBeta();
104 SelectCrossSectionTable();
105 m_isChanged = false;
106 }
107
108 return true;
109}
Garfield::TrackBichsel::GetClusterDensity
virtual double GetClusterDensity()
Definition: TrackBichsel.cc:161
Garfield::Track::GetBeta
double GetBeta() const
Definition: Track.hh:40
Garfield::RndmDirection
void RndmDirection(double &dx, double &dy, double &dz, const double length=1.)
Draw a random (isotropic) direction vector.
Definition: Random.hh:106

◆ SetDataFile()

void Garfield::TrackBichsel::SetDataFile ( const std::string &  filename)
inline

Definition at line 31 of file TrackBichsel.hh.

31{ m_datafile = filename; }
The documentation for this class was generated from the following files: