Garfield++ 3.0
A toolkit for the detailed simulation of particle detectors based on ionisation measurement in gases and semiconductors
Loading...
Searching...
No Matches
Garfield::TrackHeed Class Reference

Generate tracks using Heed++. More...

#include <TrackHeed.hh>

+ Inheritance diagram for Garfield::TrackHeed:

Public Member Functions

 TrackHeed ()
 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 &xcls, double &ycls, double &zcls, double &tcls, int &n, double &e, double &extra) override
 
bool GetCluster (double &xcls, double &ycls, double &zcls, double &tcls, int &ne, int &ni, double &e, 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
 
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).
 
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)
 
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 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)
 
- Public Member Functions inherited from Garfield::Track
 Track ()
 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.
 
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)
 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)
 
- 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
 
bool m_isElectron = false
 
std::string m_particleName = "mu-"
 
Sensorm_sensor = nullptr
 
bool m_isChanged = true
 
bool m_usePlotting = false
 
ViewDriftm_viewer = nullptr
 
bool m_debug = false
 
int m_plotId = -1
 

Detailed Description

Generate tracks using Heed++.

Definition at line 35 of file TrackHeed.hh.

Constructor & Destructor Documentation

◆ TrackHeed()

Garfield::TrackHeed::TrackHeed ( )

Constructor.

Definition at line 49 of file TrackHeed.cc.

49 : Track() {
50 m_className = "TrackHeed";
51 m_conductionElectrons.reserve(1000);
52 m_conductionIons.reserve(1000);
53
54 m_fieldMap.reset(new Heed::HeedFieldMap());
55}
std::string m_className
Definition: Track.hh:102
Track()
Constructor.
Definition: Track.cc:13
Retrieve electric and magnetic field from Sensor.
Definition: HeedFieldMap.h:15

◆ ~TrackHeed()

Garfield::TrackHeed::~TrackHeed ( )
virtual

Destructor.

Definition at line 57 of file TrackHeed.cc.

57{}

Member Function Documentation

◆ DisableDeltaElectronTransport()

void Garfield::TrackHeed::DisableDeltaElectronTransport ( )
inline

Switch simulation of delta electrons off.

Definition at line 156 of file TrackHeed.hh.

156{ m_doDeltaTransport = false; }

◆ DisableElectricField()

void Garfield::TrackHeed::DisableElectricField ( )

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

Definition at line 681 of file TrackHeed.cc.

681{ m_fieldMap->UseEfield(false); }

◆ DisableMagneticField()

void Garfield::TrackHeed::DisableMagneticField ( )

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

Definition at line 683 of file TrackHeed.cc.

683{ m_fieldMap->UseBfield(false); }

◆ EnableDeltaElectronTransport()

void Garfield::TrackHeed::EnableDeltaElectronTransport ( )
inline

Switch simulation of delta electrons on.

Definition at line 154 of file TrackHeed.hh.

154{ m_doDeltaTransport = true; }

Referenced by GarfieldPhysics::InitializePhysics().

◆ EnableElectricField()

void Garfield::TrackHeed::EnableElectricField ( )

Take the electric field into account in the stepping algorithm.

Definition at line 680 of file TrackHeed.cc.

680{ m_fieldMap->UseEfield(true); }

◆ EnableMagneticField()

void Garfield::TrackHeed::EnableMagneticField ( )

Take the magnetic field into account in the stepping algorithm.

Definition at line 682 of file TrackHeed.cc.

682{ m_fieldMap->UseBfield(true); }

◆ EnablePhotoAbsorptionCrossSectionOutput()

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

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

Definition at line 164 of file TrackHeed.hh.

164 {
165 m_usePacsOutput = on;
166 }

◆ EnablePhotonReabsorption()

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

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

Definition at line 159 of file TrackHeed.hh.

159 {
160 m_usePhotonReabsorption = on;
161 }

◆ GetCluster() [1/2]

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

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 227 of file TrackHeed.cc.

228 {
229 int ni = 0;
230 return GetCluster(xcls, ycls, zcls, tcls, n, ni, e, extra);
231}
bool GetCluster(double &xcls, double &ycls, double &zcls, double &tcls, int &n, double &e, double &extra) override
Definition: TrackHeed.cc:227

Referenced by GarfieldPhysics::DoIt(), and GetCluster().

◆ GetCluster() [2/2]

bool Garfield::TrackHeed::GetCluster ( double &  xcls,
double &  ycls,
double &  zcls,
double &  tcls,
int &  ne,
int &  ni,
double &  e,
double &  extra 
)

Definition at line 233 of file TrackHeed.cc.

235 {
236 // Initialise and reset.
237 xcls = ycls = zcls = tcls = 0.;
238 extra = 0.;
239 ne = ni = 0;
240 e = 0.;
241
242 m_deltaElectrons.clear();
243 m_conductionElectrons.clear();
244 m_conductionIons.clear();
245
246 // Make sure NewTrack has been called successfully.
247 if (!m_ready) {
248 std::cerr << m_className << "::GetCluster:\n"
249 << " Track has not been initialized. Call NewTrack first.\n";
250 return false;
251 }
252
253 if (m_particleBank.empty()) return false;
254 std::vector<Heed::gparticle*>::const_iterator end = m_particleBank.end();
255 if (m_bankIterator == end) return false;
256
257 // Look for the next cluster (i. e. virtual photon) in the list.
258 Heed::HeedPhoton* virtualPhoton = nullptr;
259 for (; m_bankIterator != end; ++m_bankIterator) {
260 // Convert the particle to a (virtual) photon.
261 virtualPhoton = dynamic_cast<Heed::HeedPhoton*>(*m_bankIterator);
262 if (!virtualPhoton) {
263 std::cerr << m_className << "::GetCluster:\n"
264 << " Particle is not a virtual photon. Program bug!\n";
265 // Try the next element.
266 continue;
267 }
268
269 // Get the location of the interaction (convert from mm to cm
270 // and shift with respect to bounding box center).
271 xcls = virtualPhoton->position().x * 0.1 + m_cX;
272 ycls = virtualPhoton->position().y * 0.1 + m_cY;
273 zcls = virtualPhoton->position().z * 0.1 + m_cZ;
274 tcls = virtualPhoton->time();
275 // Skip clusters outside the drift area or outside the active medium.
276 if (!IsInside(xcls, ycls, zcls)) continue;
277 // Add the first ion (at the position of the cluster).
278 m_conductionIons.emplace_back(
279 Heed::HeedCondElectron(Heed::point(virtualPhoton->position()), tcls));
280 ++m_bankIterator;
281 break;
282 }
283
284 // Stop if we did not find a virtual photon.
285 if (!virtualPhoton) return false;
286 // Plot the cluster, if requested.
287 if (m_usePlotting) PlotCluster(xcls, ycls, zcls);
288
289 std::vector<Heed::gparticle*> secondaries;
290 // Transport the virtual photon.
291 virtualPhoton->fly(secondaries);
292 // Get the transferred energy (convert from MeV to eV).
293 e = virtualPhoton->m_energy * 1.e6;
294
295 while (!secondaries.empty()) {
296 std::vector<Heed::gparticle*> newSecondaries;
297 // Loop over the secondaries.
298 for (auto secondary : secondaries) {
299 // Check if it is a delta electron.
300 auto delta = dynamic_cast<Heed::HeedDeltaElectron*>(secondary);
301 if (delta) {
302 extra += delta->kinetic_energy() * 1.e6;
303 const double x = delta->position().x * 0.1 + m_cX;
304 const double y = delta->position().y * 0.1 + m_cY;
305 const double z = delta->position().z * 0.1 + m_cZ;
306 if (!IsInside(x, y, z)) continue;
307 if (m_doDeltaTransport) {
308 // Transport the delta electron.
309 delta->fly(newSecondaries);
310 // Add the conduction electrons and ions to the list.
311 m_conductionElectrons.insert(m_conductionElectrons.end(),
312 delta->conduction_electrons.begin(),
313 delta->conduction_electrons.end());
314 m_conductionIons.insert(m_conductionIons.end(),
315 delta->conduction_ions.begin(),
316 delta->conduction_ions.end());
317 } else {
318 // Add the delta electron to the list, for later use.
319 deltaElectron newDeltaElectron;
320 newDeltaElectron.x = delta->position().x * 0.1 + m_cX;
321 newDeltaElectron.y = delta->position().y * 0.1 + m_cY;
322 newDeltaElectron.z = delta->position().z * 0.1 + m_cZ;
323 newDeltaElectron.t = delta->time();
324 newDeltaElectron.e = delta->kinetic_energy() * 1.e6;
325 newDeltaElectron.dx = delta->direction().x;
326 newDeltaElectron.dy = delta->direction().y;
327 newDeltaElectron.dz = delta->direction().z;
328 m_deltaElectrons.push_back(std::move(newDeltaElectron));
329 }
330 continue;
331 }
332 // Check if it is a real photon.
333 auto photon = dynamic_cast<Heed::HeedPhoton*>(secondary);
334 if (!photon) {
335 std::cerr << m_className << "::GetCluster:\n"
336 << " Particle is neither an electron nor a photon.\n";
337 }
338 extra += photon->m_energy * 1.e6;
339 const double x = photon->position().x * 0.1 + m_cX;
340 const double y = photon->position().y * 0.1 + m_cY;
341 const double z = photon->position().z * 0.1 + m_cZ;
342 if (!IsInside(x, y, z)) continue;
343 // Transport the photon.
344 if (m_usePhotonReabsorption) photon->fly(newSecondaries);
345 }
346 for (auto secondary : secondaries)
347 if (secondary) delete secondary;
348 secondaries.clear();
349 secondaries.swap(newSecondaries);
350 }
351 // Get the total number of electrons produced in this step.
352 ne = m_doDeltaTransport ? m_conductionElectrons.size()
353 : m_deltaElectrons.size();
354 ni = m_conductionIons.size();
355 return true;
356}
void PlotCluster(const double x0, const double y0, const double z0)
Definition: Track.cc:199
bool m_usePlotting
Definition: Track.hh:116
double m_energy
Photon energy [MeV].
Definition: HeedPhoton.h:36
vfloat time() const
Get the current time of the particle.
Definition: gparticle.h:176
const vec & position() const
Get the current position of the particle.
Definition: gparticle.h:174
virtual void fly(std::vector< gparticle * > &secondaries)
Transport the particle.
Definition: gparticle.h:154
double kinetic_energy() const
Get the current kinetic energy.
Definition: mparticle.h:34
Point.
Definition: vec.h:366
vfloat x
Definition: vec.h:190
vfloat z
Definition: vec.h:192
vfloat y
Definition: vec.h:191

◆ 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 195 of file TrackHeed.cc.

195 {
196 if (!m_ready) {
197 std::cerr << m_className << "::GetClusterDensity:\n"
198 << " Track has not been initialized.\n";
199 return 0.;
200 }
201
202 if (!m_transferCs) {
203 std::cerr << m_className << "::GetClusterDensity:\n"
204 << " Ionisation cross-section is not available.\n";
205 return 0.;
206 }
207
208 return m_transferCs->quanC;
209}

◆ 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 358 of file TrackHeed.cc.

360 {
361 // Make sure NewTrack has successfully been called.
362 if (!m_ready) {
363 std::cerr << m_className << "::GetElectron:\n"
364 << " Track has not been initialized. Call NewTrack first.\n";
365 return false;
366 }
367
368 if (m_doDeltaTransport) {
369 // Make sure an electron with this number exists.
370 if (i >= m_conductionElectrons.size()) {
371 std::cerr << m_className << "::GetElectron: Index out of range.\n";
372 return false;
373 }
374
375 x = m_conductionElectrons[i].x * 0.1 + m_cX;
376 y = m_conductionElectrons[i].y * 0.1 + m_cY;
377 z = m_conductionElectrons[i].z * 0.1 + m_cZ;
378 t = m_conductionElectrons[i].time;
379 e = 0.;
380 dx = dy = dz = 0.;
381
382 } else {
383 // Make sure a delta electron with this number exists.
384 if (i >= m_deltaElectrons.size()) {
385 std::cerr << m_className << "::GetElectron:\n"
386 << " Delta electron number out of range.\n";
387 return false;
388 }
389
390 x = m_deltaElectrons[i].x;
391 y = m_deltaElectrons[i].y;
392 z = m_deltaElectrons[i].z;
393 t = m_deltaElectrons[i].t;
394 e = m_deltaElectrons[i].e;
395 dx = m_deltaElectrons[i].dx;
396 dy = m_deltaElectrons[i].dy;
397 dz = m_deltaElectrons[i].dz;
398 }
399 return true;
400}

Referenced by GarfieldPhysics::DoIt().

◆ GetFanoFactor()

double Garfield::TrackHeed::GetFanoFactor ( ) const

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

Definition at line 1047 of file TrackHeed.cc.

1047{ 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 402 of file TrackHeed.cc.

403 {
404 // Make sure a "conduction" ion with this number exists.
405 if (i >= m_conductionIons.size()) {
406 std::cerr << m_className << "::GetIon: Index out of range.\n";
407 return false;
408 }
409
410 x = m_conductionIons[i].x * 0.1 + m_cX;
411 y = m_conductionIons[i].y * 0.1 + m_cY;
412 z = m_conductionIons[i].z * 0.1 + m_cZ;
413 t = m_conductionIons[i].time;
414 return true;
415}

◆ GetSteppingLimits()

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

Definition at line 145 of file TrackHeed.hh.

146 {
147 maxStep = m_maxStep;
148 radStraight = m_radStraight;
149 stepAngleStraight = m_stepAngleStraight;
150 stepAngleCurved = m_stepAngleCurved;
151 }

◆ GetStoppingPower()

double Garfield::TrackHeed::GetStoppingPower ( )
overridevirtual

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

Reimplemented from Garfield::Track.

Definition at line 211 of file TrackHeed.cc.

211 {
212 if (!m_ready) {
213 std::cerr << m_className << "::GetStoppingPower:\n"
214 << " Track has not been initialized.\n";
215 return 0.;
216 }
217
218 if (!m_transferCs) {
219 std::cerr << m_className << "::GetStoppingPower:\n"
220 << " Ionisation cross-section is not available.\n";
221 return 0.;
222 }
223
224 return m_transferCs->meanC1 * 1.e6;
225}

◆ GetW()

double Garfield::TrackHeed::GetW ( ) const

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

Definition at line 1046 of file TrackHeed.cc.

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

Referenced by GarfieldPhysics::DoIt().

◆ 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 59 of file TrackHeed.cc.

61 {
62 m_hasActiveTrack = false;
63 m_ready = false;
64
65 // Make sure the sensor has been set.
66 if (!m_sensor) {
67 std::cerr << m_className << "::NewTrack: Sensor is not defined.\n";
68 return false;
69 }
70
71 bool update = false;
72 if (!UpdateBoundingBox(update)) return false;
73
74 // Make sure the initial position is inside an ionisable medium.
75 Medium* medium = nullptr;
76 if (!m_sensor->GetMedium(x0, y0, z0, medium)) {
77 std::cerr << m_className << "::NewTrack:\n"
78 << " No medium at initial position.\n";
79 return false;
80 } else if (!medium->IsIonisable()) {
81 std::cerr << m_className << "::NewTrack:\n"
82 << " Medium at initial position is not ionisable.\n";
83 return false;
84 }
85
86 // Check if the medium has changed since the last call.
87 if (medium->GetName() != m_mediumName ||
88 fabs(medium->GetMassDensity() - m_mediumDensity) > 1.e-9) {
89 m_isChanged = true;
90 }
91
92 // If medium, particle or bounding box have changed,
93 // update the cross-sections.
94 if (m_isChanged) {
95 if (!Setup(medium)) return false;
96 m_isChanged = false;
97 m_mediumName = medium->GetName();
98 m_mediumDensity = medium->GetMassDensity();
99 }
100
101 ClearParticleBank();
102
103 m_deltaElectrons.clear();
104 m_conductionElectrons.clear();
105 m_conductionIons.clear();
106
107 // Check the direction vector.
108 double dx = dx0, dy = dy0, dz = dz0;
109 const double d = sqrt(dx * dx + dy * dy + dz * dz);
110 if (d < Small) {
111 if (m_debug) {
112 std::cout << m_className << "::NewTrack:\n"
113 << " Direction vector has zero norm.\n"
114 << " Initial direction is randomized.\n";
115 }
116 // Null vector. Sample the direction isotropically.
117 RndmDirection(dx, dy, dz);
118 } else {
119 // Normalise the direction vector.
120 dx /= d;
121 dy /= d;
122 dz /= d;
123 }
124 Heed::vec velocity(dx, dy, dz);
125 velocity = velocity * Heed::CLHEP::c_light * GetBeta();
126
127 if (m_debug) {
128 std::cout << m_className << "::NewTrack:\n Track starts at (" << x0
129 << ", " << y0 << ", " << z0 << ") at time " << t0 << "\n"
130 << " Direction: (" << dx << ", " << dy << ", " << dz << ")\n";
131 }
132
133 // Initial position (shift with respect to bounding box center and
134 // convert from cm to mm).
135 Heed::point p0((x0 - m_cX) * 10., (y0 - m_cY) * 10., (z0 - m_cZ) * 10.);
136
137 // Setup the particle.
139 if (m_particleName == "e-") {
140 particleType = &Heed::electron_def;
141 } else if (m_particleName == "e+") {
142 particleType = &Heed::positron_def;
143 } else if (m_particleName == "mu-") {
144 particleType = &Heed::muon_minus_def;
145 } else if (m_particleName == "mu+") {
146 particleType = &Heed::muon_plus_def;
147 } else if (m_particleName == "pi-") {
148 particleType = &Heed::pi_minus_meson_def;
149 } else if (m_particleName == "pi+") {
150 particleType = &Heed::pi_plus_meson_def;
151 } else if (m_particleName == "K-") {
152 particleType = &Heed::K_minus_meson_def;
153 } else if (m_particleName == "K+") {
154 particleType = &Heed::K_plus_meson_def;
155 } else if (m_particleName == "p") {
156 particleType = &Heed::proton_def;
157 } else if (m_particleName == "pbar") {
158 particleType = &Heed::anti_proton_def;
159 } else if (m_particleName == "d") {
160 particleType = &Heed::deuteron_def;
161 } else if (m_particleName == "alpha") {
162 particleType = &Heed::alpha_particle_def;
163 } else if (m_particleName == "exotic") {
164 // User defined particle
167 particleType = &Heed::user_particle_def;
168 } else {
169 // Not a predefined particle, use muon definition.
170 if (m_q > 0.) {
171 particleType = &Heed::muon_minus_def;
172 } else {
173 particleType = &Heed::muon_plus_def;
174 }
175 }
176
177 Heed::HeedParticle particle(m_chamber.get(), p0, velocity, t0, particleType,
178 m_fieldMap.get());
179 // Set the step limits.
180 particle.set_step_limits(m_maxStep * Heed::CLHEP::cm,
181 m_radStraight * Heed::CLHEP::cm,
182 m_stepAngleStraight * Heed::CLHEP::rad,
183 m_stepAngleCurved * Heed::CLHEP::rad);
184 // Transport the particle.
185 particle.fly(m_particleBank);
186 m_bankIterator = m_particleBank.begin();
187 m_hasActiveTrack = true;
188 m_ready = true;
189
190 // Plot the new track.
191 if (m_usePlotting) PlotNewTrack(x0, y0, z0);
192 return true;
193}
bool GetMedium(const double x, const double y, const double z, Medium *&medium)
Get the medium at (x, y, z).
Definition: Sensor.cc:166
Sensor * m_sensor
Definition: Track.hh:112
bool m_debug
Definition: Track.hh:119
double GetBeta() const
Return the speed ( ) of the projectile.
Definition: Track.hh:54
bool m_isChanged
Definition: Track.hh:114
double m_q
Definition: Track.hh:104
std::string m_particleName
Definition: Track.hh:110
double m_mass
Definition: Track.hh:106
void PlotNewTrack(const double x0, const double y0, const double z0)
Definition: Track.cc:193
void set_mass(const double m)
void set_charge(const double z)
Definition: vec.h:177
void RndmDirection(double &dx, double &dy, double &dz, const double length=1.)
Draw a random (isotropic) direction vector.
Definition: Random.hh:107
particle_def pi_minus_meson_def("pi_minus_meson", "pi-", 139.56755 *MeV/c_squared, -eplus, 0, 0, 0.0, spin_def(1.0, -1.0))
Definition: particle_def.h:114
particle_def pi_plus_meson_def("pi_plus_meson", "pi+", 139.56755 *MeV/c_squared, eplus, 0, 0, 0.0, spin_def(1.0, 1.0))
Definition: particle_def.h:112
particle_def muon_minus_def("muon_minus", "mu-", 105.658367 *MeV/c_squared, electron_charge, 1, 0, 0.5, spin_def(0.0, 0.0))
Definition: particle_def.h:101
particle_def alpha_particle_def("alpha_particle", "alpha", 3727.417 *MeV/c_squared, 2 *eplus, 0, 4, 0.0, spin_def(0.0, 0.0))
Definition: particle_def.h:120
particle_def anti_proton_def("", "p-", proton_def)
Definition: particle_def.h:104
particle_def deuteron_def("deuteron", "dtr", 1875.613 *MeV/c_squared, eplus, 0, 2, 0.0, spin_def(0.0, 0.0))
Definition: particle_def.h:119
particle_def proton_def("proton", "p+", proton_mass_c2/c_squared, eplus, 0, 1, 0.5, spin_def(0.5, 0.5))
Definition: particle_def.h:103
particle_def K_minus_meson_def("K_minus_meson_def", "K-", K_plus_meson_def)
Definition: particle_def.h:117
particle_def K_plus_meson_def("K_plus_meson_def", "K+", 493.677 *MeV/c_squared, 1, 0, 0, 0.0, spin_def(0.5, -0.5))
Definition: particle_def.h:116
particle_def user_particle_def("user_particle", "X", 139.56755 *MeV/c_squared, eplus, 0, 0, 0.0, spin_def(0.0, 0.0))
Definition: particle_def.h:123
DoubleAc fabs(const DoubleAc &f)
Definition: DoubleAc.h:615
particle_def muon_plus_def("muon_plus", "mu+", muon_minus_def)
Definition: particle_def.h:102
particle_def positron_def("positron", "e+", electron_def)
Definition: particle_def.h:100
particle_def electron_def("electron", "e-", electron_mass_c2/c_squared, electron_charge, 1, 0, 0.5, spin_def(0.0, 0.0))
Definition: particle_def.h:99
DoubleAc sqrt(const DoubleAc &f)
Definition: DoubleAc.cpp:314

Referenced by GarfieldPhysics::DoIt().

◆ 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 685 of file TrackHeed.cc.

686 {
687 if (fabs(e1 - e0) < Small) {
688 std::cerr << m_className << "::SetEnergyMesh:\n"
689 << " Invalid energy range:\n"
690 << " " << e0 << " < E [eV] < " << e1 << "\n";
691 return;
692 }
693
694 if (nsteps <= 0) {
695 std::cerr << m_className << "::SetEnergyMesh:\n"
696 << " Number of intervals must be > 0.\n";
697 return;
698 }
699
700 m_emin = std::min(e0, e1);
701 m_emax = std::max(e0, e1);
702 m_emin *= 1.e-6;
703 m_emax *= 1.e-6;
704 m_nEnergyIntervals = nsteps;
705}

◆ 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 707 of file TrackHeed.cc.

707 {
708 if (fabs(z) < Small) {
709 std::cerr << m_className << "::SetParticleUser:\n"
710 << " Particle cannot have zero charge.\n";
711 return;
712 }
713 if (m < Small) {
714 std::cerr << m_className << "::SetParticleUser:\n"
715 << " Particle mass must be greater than zero.\n";
716 }
717 m_q = z;
718 m_mass = m;
719 m_isElectron = false;
720 m_spin = 0;
721 m_particleName = "exotic";
722}
bool m_isElectron
Definition: Track.hh:109

◆ 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 137 of file TrackHeed.hh.

139 {
140 m_maxStep = maxStep;
141 m_radStraight = radStraight;
142 m_stepAngleStraight = stepAngleStraight;
143 m_stepAngleCurved = stepAngleCurved;
144 }

◆ 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 417 of file TrackHeed.cc.

421 {
422 int ni = 0;
423 return TransportDeltaElectron(x0, y0, z0, t0, e0, dx0, dy0, dz0, nel, ni);
424}
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:426

◆ 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 426 of file TrackHeed.cc.

430 {
431 nel = 0;
432 ni = 0;
433
434 // Check if delta electron transport was disabled.
435 if (!m_doDeltaTransport) {
436 std::cerr << m_className << "::TransportDeltaElectron:\n"
437 << " Delta electron transport has been switched off.\n";
438 return;
439 }
440
441 // Make sure the sensor has been set.
442 if (!m_sensor) {
443 std::cerr << m_className << "::TransportDeltaElectron:\n"
444 << " Sensor is not defined.\n";
445 m_ready = false;
446 return;
447 }
448
449 bool update = false;
450 if (!UpdateBoundingBox(update)) return;
451
452 // Make sure the initial position is inside an ionisable medium.
453 Medium* medium = nullptr;
454 if (!m_sensor->GetMedium(x0, y0, z0, medium)) {
455 std::cerr << m_className << "::TransportDeltaElectron:\n"
456 << " No medium at initial position.\n";
457 return;
458 } else if (!medium->IsIonisable()) {
459 std::cerr << "TrackHeed:TransportDeltaElectron:\n"
460 << " Medium at initial position is not ionisable.\n";
461 m_ready = false;
462 return;
463 }
464
465 // Check if the medium has changed since the last call.
466 if (medium->GetName() != m_mediumName ||
467 fabs(medium->GetMassDensity() - m_mediumDensity) > 1.e-9) {
468 m_isChanged = true;
469 update = true;
470 m_ready = false;
471 m_hasActiveTrack = false;
472 }
473
474 // If medium or bounding box have changed, update the "chamber".
475 if (update) {
476 if (!Setup(medium)) return;
477 m_ready = true;
478 m_mediumName = medium->GetName();
479 m_mediumDensity = medium->GetMassDensity();
480 }
481
482 m_deltaElectrons.clear();
483 m_conductionElectrons.clear();
484 m_conductionIons.clear();
485
486 // Initial position (shift with respect to bounding box center and
487 // convert from cm to mm).
488 Heed::point p0((x0 - m_cX) * 10., (y0 - m_cY) * 10., (z0 - m_cZ) * 10.);
489
490 // Make sure the kinetic energy is positive.
491 if (e0 <= 0.) {
492 // Just create a conduction electron on the spot.
493 m_conductionElectrons.emplace_back(Heed::HeedCondElectron(p0, t0));
494 nel = 1;
495 return;
496 }
497
498 // Check the direction vector.
499 double dx = dx0, dy = dy0, dz = dz0;
500 const double d = sqrt(dx * dx + dy * dy + dz * dz);
501 if (d <= 0.) {
502 // Null vector. Sample the direction isotropically.
503 RndmDirection(dx, dy, dz);
504 } else {
505 // Normalise the direction vector.
506 dx /= d;
507 dy /= d;
508 dz /= d;
509 }
510 Heed::vec velocity(dx, dy, dz);
511
512 // Calculate the speed for the given kinetic energy.
513 const double gamma = 1. + e0 / ElectronMass;
514 const double beta = sqrt(1. - 1. / (gamma * gamma));
515 double speed = Heed::CLHEP::c_light * beta;
516 velocity = velocity * speed;
517
518 // Transport the electron.
519 std::vector<Heed::gparticle*> secondaries;
520 Heed::HeedDeltaElectron delta(m_chamber.get(), p0, velocity, t0, 0,
521 m_fieldMap.get());
522 delta.fly(secondaries);
523 ClearBank(secondaries);
524
525 m_conductionElectrons.swap(delta.conduction_electrons);
526 m_conductionIons.swap(delta.conduction_ions);
527 nel = m_conductionElectrons.size();
528 ni = m_conductionIons.size();
529}

Referenced by GarfieldPhysics::DoIt(), and TransportDeltaElectron().

◆ TransportPhoton() [1/2]

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 531 of file TrackHeed.cc.

534 {
535 int ni = 0;
536 TransportPhoton(x0, y0, z0, t0, e0, dx0, dy0, dz0, nel, ni);
537}
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)
Definition: TrackHeed.cc:539

◆ TransportPhoton() [2/2]

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
ne,ninumber of electrons/ions produced by the photon

Definition at line 539 of file TrackHeed.cc.

543 {
544 nel = 0;
545 ni = 0;
546
547 // Make sure the energy is positive.
548 if (e0 <= 0.) {
549 std::cerr << m_className << "::TransportPhoton:\n"
550 << " Photon energy must be positive.\n";
551 return;
552 }
553
554 // Make sure the sensor has been set.
555 if (!m_sensor) {
556 std::cerr << m_className << "::TransportPhoton: Sensor is not defined.\n";
557 m_ready = false;
558 return;
559 }
560
561 bool update = false;
562 if (!UpdateBoundingBox(update)) return;
563
564 // Make sure the initial position is inside an ionisable medium.
565 Medium* medium = nullptr;
566 if (!m_sensor->GetMedium(x0, y0, z0, medium)) {
567 std::cerr << m_className << "::TransportPhoton:\n"
568 << " No medium at initial position.\n";
569 return;
570 } else if (!medium->IsIonisable()) {
571 std::cerr << "TrackHeed:TransportPhoton:\n"
572 << " Medium at initial position is not ionisable.\n";
573 m_ready = false;
574 return;
575 }
576
577 // Check if the medium has changed since the last call.
578 if (medium->GetName() != m_mediumName ||
579 fabs(medium->GetMassDensity() - m_mediumDensity) > 1.e-9) {
580 m_isChanged = true;
581 update = true;
582 m_ready = false;
583 }
584
585 // If medium or bounding box have changed, update the "chamber".
586 if (update) {
587 if (!Setup(medium)) return;
588 m_ready = true;
589 m_mediumName = medium->GetName();
590 m_mediumDensity = medium->GetMassDensity();
591 }
592
593 // Delete the particle bank.
594 // Clusters from the current track will be lost.
595 m_hasActiveTrack = false;
596 ClearParticleBank();
597 m_deltaElectrons.clear();
598 m_conductionElectrons.clear();
599 m_conductionIons.clear();
600
601 // Check the direction vector.
602 double dx = dx0, dy = dy0, dz = dz0;
603 const double d = sqrt(dx * dx + dy * dy + dz * dz);
604 if (d <= 0.) {
605 // Null vector. Sample the direction isotropically.
606 RndmDirection(dx, dy, dz);
607 } else {
608 // Normalise the direction vector.
609 dx /= d;
610 dy /= d;
611 dz /= d;
612 }
613 Heed::vec velocity(dx, dy, dz);
614 velocity = velocity * Heed::CLHEP::c_light;
615
616 // Initial position (shift with respect to bounding box center and
617 // convert from cm to mm).
618 Heed::point p0((x0 - m_cX) * 10., (y0 - m_cY) * 10., (z0 - m_cZ) * 10.);
619
620 // Create and transport the photon.
621 Heed::HeedPhoton photon(m_chamber.get(), p0, velocity, t0, 0, e0 * 1.e-6,
622 m_fieldMap.get());
623 std::vector<Heed::gparticle*> secondaries;
624 photon.fly(secondaries);
625
626 while (!secondaries.empty()) {
627 std::vector<Heed::gparticle*> newSecondaries;
628 // Loop over the particle bank and look for daughter particles.
629 std::vector<Heed::gparticle*>::iterator it;
630 for (it = secondaries.begin(); it != secondaries.end(); ++it) {
631 // Check if it is a delta electron.
632 auto delta = dynamic_cast<Heed::HeedDeltaElectron*>(*it);
633 if (delta) {
634 if (m_doDeltaTransport) {
635 // Transport the delta electron.
636 delta->fly(newSecondaries);
637 // Add the conduction electrons to the list.
638 m_conductionElectrons.insert(m_conductionElectrons.end(),
639 delta->conduction_electrons.begin(),
640 delta->conduction_electrons.end());
641 m_conductionIons.insert(m_conductionIons.end(),
642 delta->conduction_ions.begin(),
643 delta->conduction_ions.end());
644 } else {
645 // Add the delta electron to the list, for later use.
646 deltaElectron newDeltaElectron;
647 newDeltaElectron.x = delta->position().x * 0.1 + m_cX;
648 newDeltaElectron.y = delta->position().y * 0.1 + m_cY;
649 newDeltaElectron.z = delta->position().z * 0.1 + m_cZ;
650 newDeltaElectron.t = delta->time();
651 newDeltaElectron.e = delta->kinetic_energy() * 1.e6;
652 newDeltaElectron.dx = delta->direction().x;
653 newDeltaElectron.dy = delta->direction().y;
654 newDeltaElectron.dz = delta->direction().z;
655 m_deltaElectrons.push_back(std::move(newDeltaElectron));
656 }
657 continue;
658 }
659 // Check if it is a fluorescence photon.
660 auto fluorescencePhoton = dynamic_cast<Heed::HeedPhoton*>(*it);
661 if (!fluorescencePhoton) {
662 std::cerr << m_className << "::TransportPhoton:\n"
663 << " Unknown secondary particle.\n";
664 ClearBank(secondaries);
665 ClearBank(newSecondaries);
666 return;
667 }
668 fluorescencePhoton->fly(newSecondaries);
669 }
670 secondaries.swap(newSecondaries);
671 ClearBank(newSecondaries);
672 }
673 ClearBank(secondaries);
674 // Get the total number of electrons produced in this step.
675 nel = m_doDeltaTransport ? m_conductionElectrons.size()
676 : m_deltaElectrons.size();
677 ni = m_conductionIons.size();
678}

Referenced by GarfieldPhysics::DoIt(), and TransportPhoton().


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